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I.  INTRODUCTION 


The  objective  of  this  work  is  the  investigation  of  vortical  flows  over  three- 
dimensional  bodies  at  high  incidence  utilizing  numerical  methods.  The  advantage  of 
numerical  simulations  compared  with  experiment  is  that  they  allow  simultaneous  obser¬ 
vation  of  all  flow  quantities  of  interest  for  the  entire  flow  field.  The  disadvantage  of 
numerical  techniques  is  the  accuracy  limitations  for  simulations  of  flow  fields  over 
complex  realistic  configurations, .even  with  the  most  efficient  numerical  schemes  and  fast 
computers.  High  Reynolds  number  turbulent  flows  of  engineering  interest  can  be  fully 
simulated  when  all  relevant  scales  are  resolved.  Resolution  of  all  scales  for  complex  high 
reynolds  number  flows  over  realistic  configurations  is  beyond  the  capabilities  of  the 
present  and  next  generation  supercomputers.  Common  practice  for  the  simulation  of 
engineering  flows  is  the  use  of  various  turbulence  models  to  approximate  the  effect  of 
the  small  scales  which  cannot  be  resolved.  Error  sources  in  numerical  simulations  are 
related  to  the  discretization  process,  the  order  of  accuracy  of  the  numerical  scheme  and 
the  turbulence  modeling  that  is  used. 

Nevertheless.  Computational  Fluid  Dynamics  (CFD)  allows  investigation  of  various 
fluid  How  phenomena  that  in  the  past  was  possible  only  in  wind  tunnels,  water  tunnels 
or  actual  flight  testing.  The  advantage  of  being  able  to  accurately  capture  the  flow 
characteristics  over  complex  configurations  or  even  complete  aircraft  without  endan¬ 
gering  life,  i.e.  preliminary  flight  testing,  is  readily  apparent.  Numerical  solutions  also 
enable  to  investigate  and  visualize  the  flow  field  characteristics  from  any  viewpoint  or 
in  as  much  detail  as  desired.  With  the  ever  increasing  speed  cost  ratio  of  today's  com¬ 
puters.  CFD  techniques  will  be  playing  a  more  significant  role  facilitating  aerodynamic 
research  and  supplementing  experimental  investigations.  Even  though  CFD  and 
Navier-Stokes  methods  are  not  a  new  research  tool,  new  and  more  efficient  numerical 
techniques  are  evolving,  while  at  the  same  time  computers  are  becoming  faster.  Nu¬ 
merical  prediction  of  steady  flows  over  complete  aircraft  and  comparison  with  flight  data 
is  already  undent  ay  [Ref.  1J.  In  the  near  future  CFD  is  expected  to  play  a  more  active 
role  in  fluid  dynamic  research  enabling  simulation  of  complex  unsteady  flow  regimes. 

In  the  past  panel  methods  and  vortex  lattice  methods  were  used  in  the  analysis  of 
flows  [Ref.  2].  These  methods  were  insufficient  for  a  detailed  analysis  of  complex  flows 
such  as  vortical  flows  over  bodies  at  incidence.  The  limitations  of  these  methods  are  due 
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to  the  potential  flow  assumption  which  is  valid  only  for  inviscid  and  irrotational  flow. 
Viscous  effects  close  to  the  surface  for  attached  or  mildly  separated  flow  are  obtained 
using  Boundary  Layer  methods  [Ref.  3].  The  rotational  compressible  flow  regime  at  high 
Reynolds  numbers  was  investigated  with  the  Euler  equations.  Viscous  effects  become 
more  important  for  flows  at  high  angles  of  attack;  therefore,  the  solution  of  the  Navier- 
Stokes  equations  is  required  for  this  flow  regime. 

In  Chapter  2  the  theoretical  development  of  the  compressible  Navier-Stokes 
equations  will  be  discussed.  A  finite  difference  algorithm  used  for  the  numerical  solution 
of  these  equations  will  be  presented.  The  numerical  solution  is  performed  on  the  finite 
number  of  points  obtained  after  discretization  of  the  flow  domain.  The  procedure  which 
yields  this  finite  collection  of  points  in  the  solution  domain  is  known  as  grid  generation. 
The  quality  of  the  solution  depends  directly  on  the  smoothness  of  the  grid  and  its  ability 
to  accurately  represent  flow  gradients.  Therefore,  grid  generation  is  an  important  part 
of  the  numerical  solution.  However,  the  numerical  solution  of  the  governing  equations 
is  not  the  main  objective  of  this  research.  The  grid  generation  part,  which  is  a  necessary 
stage  before  starting  the  numerical  solution  will  be  covered  in  full  detail.  Numerical 
solutions  depend  on  the  representation  of  the  flow  field  by  an  orderly,  finite  collection 
of  points.  The  process  of  obtaining  three-dimensional  grids  involves  first  definition  of 
an  inner  boundary,  commonly  known  as  the  surface  grid,  before  the  subsequent  gener¬ 
ation  of  the  field  grid  can  begin. 

The  methods  available  for  both  the  surface  and  field  grid  generation  will  be  covered 
in  Chapters  3  and  4.  respectively.  Developments  in  the  area  of  grid  generation  have 
provided  a  key  to  eliminate  the  problem  of  boundary  shape  definition  [Ref.  4j.  Finite 
difference  grids  can  also  be  used  to  construct  meshes  that  are  suitable  in  finite  element 
methods.  The  specific  numerical  method  utilized  in  this  research  is  the  finite  difference 
method.  The  finite  difference  method  is  one  of  the  oldest  numerical  methods  that  can 
be  utilized  to  obtain  numerical  solutions  to  differential  equations.  The  application  of 
this  method  is  based  on  a  Taylor  series  expansion  and  the  definition  of  the  derivative: 
most  likely  first  developed  by  Euler  in  176S  [Ref.  5;  p.  167],  The  algorithm  used  for  the 
numerical  integration  utilizes  a  partially  flux-split  numerical  scheme  with  central  differ¬ 
encing  in  the  other  two  directions  [Ref.  6j. 

The  methods  described  above  will  be  applied  to  a  double-delta  wing  that  has  a  strake 
with  a  sweep  angle  of  76°  and  a  delta  wing  section  with  a  sweep  angle  of  40° .  Particular 
emphasis  will  be  placed  on  the  investigation  of  the  vortical  flow  field  at  moderate  to  high 
angles  of  attack.  Separated  flow  along  the  strake's  leading  edge  forms  free  shear  layers 
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which  roll  up  to  form  vortex  cores.  This  primary  strake  vortex  generates  an  additional 
non-linear  lift  called  vortex  induced  lift.  A  primary  wing  vortex  also  develops  from  the 
leading  edge  of  the  40°  swept  delta  wing. 

The  mutual  interaction  of  the  strake  and  wing  vortices  and  their  interaction  with  the 
surface  is  an  active  area  of  current  research.  Investigation  and  prediction  of  the  vortex 
breakdown  that  appears  at  higher  angles  of  attack  is  also  of  high  interest.  The  devel¬ 
opment  of  the  leading  edge  vortex  as  well  as  breakdown  are  important  phenomena  that 
need  to  be  fully  understood.  Various  angles  of  attack,  a  =  10.0° ,  19.0° ,  22.4°  are  in¬ 
vestigated  and  compared  with  available  experimental  data.  Understanding  the  leeward- 
side  flow  structure  as  well  as  breakdown  are  important  phenomena  that  affect 
significantly  today's  tactical  and  fighter  aircraft  effectiveness. 

Vortex  breakdown  is  a  transition  of  the  vortex  core  from  a  jet-like  flow  to  a  wake- 
like  flow.  Both  swirl  angle  and  adverse  pressure  gradient  along  the  axial  direction  con¬ 
tribute  to  the  breakdown  of  the  vortex.  Peckham  and  Atkinson  first  identified  vortex 
breakdown  when  analyzing  delta  wings  at  high  angles  of  attack  [Ref.  7],  Research  on 
vortex  breakdown  was  continued  by  Elle,  Lambourne  and  Bryer,  Harvey.  Pritchard, 
Sarpkaya.  Hummel,  Faler  and  Leibovitch,  Payne  and  Nelson  [Ref.  8,9,10,11,  12 
.13,14.15,16].  Studies  then  naturally  progressed  to  more  complicated  bodies  such  as  the 
double-delta  wing  where  Brennenstuhl  tested  several  wings  in  a  low  speed  wind  tunnel 
and  a  water  tunnel  [Ref.  17],  The  present  study  will  attempt  a  comparison  of  the  com¬ 
putational  solution  with  the  data  obtained  from  wind  tunnel  testing  done  by 
Cunningham  and  Boer  [Ref.  IS].  This  comparison  along  with  discussions  of  the  results 
that  were  developed  during  this  research  will  be  covered  in  Chapter  5.  The  closing 
chapter  summarizes  the  conclusions  and  presents  recommendations  for  further  research. 


II.  THEORETICAL  APPROACH 


The  main  objective  of  this  work  is  the  investigation  of  different  techniques  for  the 
grid  generation  over  complex  three-dimensional  bodies,  and  numerical  flow  visualization 
of  the  computed  flowfields  over  bodies  at  high  angles  of  attack.  The  flow  field  is  ob¬ 
tained  by  the  numerical  solution  of  the  Navier-Stokes  equations.  Fluid  flow  in  the 
continuum  (low  regime  includes  most  of  the  physical  flows  and  is  governed  by  the 
Navier-Stokes  equations.  The  derivation  of  the  Navier-Stokes  equations  is  well  known 
[Ref.  3:  pp.  47-66).  Solutions  of  these  equations  are  of  interest  in  basic  fluid  mechanics 
research  and  for  engineering  applications.  The  solution  of  the  Navier-Stokes  equations 
is  quite  difficult  due  to  their  nonlinearity.  Analytical  closed  form  solutions  of  the 
Navier-Stokes  equations  can  be  obtained  for  only  a  few  flow  situations  of  simple  ge¬ 
ometrical  configurations  and  boundary  conditions.  Simplified  forms  of  the  Navier- 
Stokes  equations,  such  as  the  boundary  layer  equations,  can  give  satisfactory  answers 
for  many  flows  of  practical  interest.  The  inviscid  form  of  the  Navier-Stokes  equations, 
commonly  known  as  the  Euler  equations,  can  provide  solutions  for  flows  away  from 
solid  boundaries.  However,  complex  flows  such  as  vortical  separated  flows  require  the 
solution  of  the  full  Navier-Stokes  equations,  which  can  only  be  obtained  by  utilizing 
numerical  techniques.  The  derivation  of  the  Navier-Stokes  equations  is  outlined  in  the 
following  paragraphs. 

A.  GOVERNING  EQUATIONS 
1.  The  Continuity  Equation 

For  the  derivation  of  the  Navier-Stokes  equations  the  fluid  medium  is  consid¬ 
ered  as  an  isotropic,  homogeneous,  compressible  and  viscous  Newtonian  fluid.  The 
continuity  equation  is  a  manifestation  of  the  fact  that  mass  can  neither  be  created  nor 
destroyed.  The  continuity  equation  states  that  the  time  variation  of  density  within  a 
control  volume  plus  the  mass  entering  and  leaving  the  control  volume  is  equal  to  zero. 
The  differential  form  of  the  continuity  equation  for  a  compressible  fluid  and  noivsteady 
flow  can  be  written  as; 


4r-  + V.(p?r)  =  0. 

ci  v  ' 

For  low  speeds  the  density  variation  is  small,  therefore; 


(1) 
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cp 

~T~  =  0 
Cl 


and, 


?  .  (p  V)  =  pv  .  f. 


Hence,  for  incompressible  flow,  the  continuity  equation  can  be  written  as; 

V  .  V  =  0. 

2.  Derivation  of  the  Navier-Stokes  Equations 

For  a  compressible  fluid,  all  primitive  variables,  density  (p),  pressure  (p)  and  velocity 
(F)  are  functions  of  space  and  time.  Newton's  Second  Law  states  that  the  summation 
of  all  forces  must  be  equal  to  the  mass  times  the  acceleration. 

I F  =  Ma  (2) 


Considering  an  infinitesimally  small  fluid  particle  or  control  volume  moving  in  a 
Cartesian  Coordinate  System,  the  right-hand  side  of  equation  (1)  can  be  rewritten  as; 


*/««  ^{pV)dxdydz 


(3) 


where  ~~  is  the  material  derivative. 

D/ 


D  fC  D  p  ,  D 

dT^’^d?1  +u~dTp 


and  F  is  the  velocity  vector,  which  for  a  Cartesian  Coordinate  system  is, 


l '  =  ui  +  vj  +  wk  (4) 

here  u,  v.  w  are  the  velocity  components  along  the  coordinate  axes.  The  external  forces 
normally  consist  of  the  gravitational  forces  and  the  forces  acting  on  the  boundaries  of 
the  control  volume,  namely  pressure  and  friction.  All  other  body  forces,  such  as 
electromagnetic  forces  will  be  ignored.  For  simplicity,  the  momentum  equation  only  for 
the  x-direction  will  be  derived,  while  the  derivation  is  analogous  for  the  y  and  z- 
directions.  The  x-componcnt  of  equation  (2)  is; 
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M  «x  =  -p7  (pu)dxdydz.  (5) 

Figure  (1)  shows  the  normal  and  shear  stresses  on  an  infinitesimal  control  volume. 
Summation  of  the  forces  in  the  x-direction  yields; 

C  T  fit 

Surface  Forces  =  oyiydi  +  (ryx  +  dy)dxdz  +  (t2x  H — r~  dz)dxdy 


cnx 

~  K  +  dxjdyd:  -  t2Xdxdy  -  ryxdxdz 


which  reduces  to. 


,  COx  ?T.xv  c r.,  , 

Surface  Forces  —  ( - r~“  H — r—  H — r-  )dxdvdz 

cx  cv  c: 


Among  the  body  forces  only  the  gravity  will  be  considered.  Therefore,  if  f  is  the  x- 
component  of  the  gravity  force  then; 

U 'eight  =  ffx  y.z)p{xy,z  )dxdydz .  (7) 

The  sum  of  equations  (6)  and  (7)  are  the  external  forces  which  are  equal  to  the  acceler¬ 
ation  as  stated  by  equation  (5).  After  cancellation  of  the  common  term  of  volume 
(dxdydz).  the  following  force  balance  for  the  x-direction  is  obtained. 


n  ,  CC 7.  CTXV  c T.,  x 

fT"  ( PU )  —  pfx  +  (  —  “~7 -  -I - “  +  — 7.  ) 

Dr  v  cx  cv  cz  ' 


The  next  step  is  to  express  the  stresses  in  terms  of  the  primitive  variables,  i.e.,  velocities 
and  pressure.  First,  the  static  pressure  is  defined  as  the  mean  of  the  normal  stresses. 

P  “  y  (°x  +  ay  + 

This  equation  can  be  algebraically  rewritten  as; 


°x=P  +  —  ~  °y  - 
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Figure  I.  Normal  and  Shear  Stress  due  to  Friction 

In  equation  (9)  the  left-hand  side  is  the  normal  stress  at  a  point  in  the  fluid.  The  first 
right-hand  side  term  is  the  static  pressure  and  the  second  right-hand  side  term  is  the 
deviation  of  the  normal  stress  from  the  pressure  due  to  viscous  forces.  Next  a  relation¬ 
ship  between  stress  and  rate  of  strain  must  be  found.  Isotropy  implies  that  this  relation 
between  the  components  of  stress  and  rate  of  strain  is  the  same  for  every  direction.  The 
Newtonian  fluid  assumption  means  that  this  relationship  is  also  linear.  Referring  to 
Figure  (2)  where  o,  and  ot  arc  resolved  into  diagonal  components  and  equating  the 
forces,  the  following  force  balance  is  obtained. 

t'joK  2  )  +  ox{  -4=-  )  -  oy{  )  =  0 
\  *  \  - 

This  equation  can  be  rewritten  as: 

-t',v  =  tK-0'  0<» 
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Figure  2.  Normal  and  Shear  Stress 


A  similar  equation  can  be  derived  for  the  xz-plane. 


T  2X  ~  2 


(ID 


Substitution  of  equations  (10)  and  (11)  into  equation  (9)  results  in; 

°X  =  P  ^  (T  2X  ~  T  JC>’)' 


(12) 


The  deformation  of  the  initial  shape  of  the  fluid  element  (ABCD)  to  (A  BC  D)  due  to 
stresses  in  the  y-direction  is  shown  in  Figure  (3).  The  same  figure  also  shows  that  the 
length  of  OA'  is  as  follows; 

(O.l')  =  Length  =  (u  +  +  H.O.T.)kt. 

The  length  change  due  to  stress  is; 
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Figure  3.  Rate  of  Strain 


{AA')  =  A  Length  s: 


cu 


cx  -2 


The  strain  is  the  change  in  length  divided  by  the  original  length  which  is  a/N  2  .  The 
strain  rate  is  obtained  by  dividing  this  length  by  A i  .  The  same  procedure  can  be  re¬ 
peated  for  the  y-direction  to  obtain  the  resultant  rate  of  strain  on  a  45  degree  plane  due 
to  a,  .  The  shear  stress  on  the  45°  plane  due  to  a,  and  oy  is  given  by  equation  (13).  A 
similar  procedure  for  the  zx-plane  yields  the  analogous  shear  stress  which  is  shown  in 
equation  ( 14). 


,  (  ru  <  v  \ 


cr 


cx 


cv 


(13) 


(14) 


2X 


-*( 


civ 


du 


cz 


cx 


In  the  above  equations,  the  proportionality  constant  p,  is  defined  as  the  coefficient  of 
viscosity.  Substituting  equation  (13)  and  (14)  back  into  equation  (12),  and  using  Stokes' 
Hypothesis; 


3/  +  2p  =  0 


yields  equation  (15)  [Ref.  3:  pp.  60-61]. 


Cll  Z  /  IU  ,  Cl  .  CIV  \\ 

-p-n{  2- - 7-(  —  +  —  +  -r- )) 


CU 


ci- 


ow 


CX 


cx 


cy 


cz 


(15) 


In  this  equation  the  first  term  in  parenthesis  is  the  linear  strain  rate  and  the  second  term 
is  the  volumetric  strain  rate.  To  complete  the  derivation,  the  terms  r,,  and  t„  in 
equation  (S  )  will  be  expressed  in  terms  of  the  velocity  components.  From  Figure  (4  )  the 
rate  of  strain  can  be  obtained  [Ref.  19:  p.  93]. 

The  rate  of  strain  on  this  element  is  Assuming  that  the  variation  of  the 

rate  of  strain  (•/)  is  small,  the  following  expressions  can  be  written. 


A;- 


Ay 


cu 

CY 


AvA  r 


Ay 


+ 


cv 

c.r 


AxAl 


Ax 


a  (  ciJ  ,  <~r  a  . 

Ay  =  ( —  +  —  )A/ 

cy  cx 


Ay 

Taking  the  limit  of  as  At  -*  0  ,  the  rate  of  strain  is  given  by; 


dy_ 

di 


cu  cv 
x*  cy  cx 


(16) 


Due  to  isotropy.  is  equal  to  Tyt  .  Analogous  procedures  used  to  derive  equation  (16) 
can  be  repeated  for  the  yz-plane  and  zx-plane,  respectively,  so  that  for  a  Newtonian 
Fluid  equations  (17)  and  (IS)  can  be  written- 
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Figure  4.  Rate  of  Strain 


=  ^  +  SL) 


'■xy 


cv  cx 


(  c'v  .  Cl/  ^ 

T~  =  ^17 +  7T} 


(17) 

(is) 


Finally,  substituting  equations  (15).  (17)  and  (18)  back  into  equation  (8),  the  momentum 
equation  for  the  x-direction  is  obtained. 


Similarly,  the  momentum  equations  for  the  y-direction  and  z-direction  can  be  derived. 
These  equations  are; 


and, 


(20) 


(21) 


The  unknowns  in  these  last  three  equations  are  the  primary  variables;  the  density,  the 
velocities  and  the  pressure,  (p,  u.  v,  w,  />)  .  The  momentum  equations  along  with  the 
continuity  equation  constitute  the  Navier-Stokes  equations  in  the  primitive  variable 
formulation.  The  continuity  equation  for  a  Cartesian  coordinate  system  is  restated; 

I±  +  ££1L  +  S£L  +  I^L,0.  (22) 

ci  cx  cy  cz 


Here  the  pressure  is  related  to  the  density  through  the  equation  of  state. 

p  —  pRT  =0  (23) 

For  an  isothermal  process  and  incompressible  flow,  equation  (19)  through  equation  (23) 
would  be  sufficient,  but  when  temperature  variations  depend  on  density  and  pressure, 
the  energy  equation  is  also  required.  This  is  always  the  case  for  compressible  flow  where 
density  depends  on  pressure  and  temperature.  The  energy  equation  expresses  the  bal¬ 
ance  between  heat  and  mechanical  energy.  The  variation  of  viscosity  due  to  temperature 
variation  may  be  obtained  by  an  empirical  viscosity  law.  The  final  result  is  a  system  of 
five  partial  differential  equations  with  five  unknowns;  u ,  v,  w,  p  ,  and  p  .[Ref.  3,  19  ] 

3.  Derivation  of  the  Energy  Equation 

It  is  well  known  that  energy  can  be  neither  created  nor  destroyed  but  it  can  only 
change  in  form.  Therefore  an  energy  balance  exists  for  a  fluid  element  in  motion.  This 
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energy  balance  is  obtained  through  certain  mechanisms  which  for  a  compressible  fluid 
are  determined  by  changes  in  heat  content,  total  energy  and  mechanical  work.  Changes 
in  heat  content  can  be  due  to  convection,  conduction,  friction  and/or  radiation.  For  the 
following  derivation  radiation  is  neglected  because  its  effect  is  small  at  moderate  tem¬ 
peratures.  The  energy  balance  for  a  control  volume  is  expressed  by  the  first  law  of 
thermodynamics. 


dQ_dE  dlV 

di  dt  dt 


(24) 


First  the  variation  of  mechanical  work  or  the  contribution  to  work  done  by  the  external 
forces  acting  along  the  x-direction  is  derived.  Again  referring  to  Figure  (1),  the  con¬ 
tribution  to  work  done  by  each  of  the  stress  components  is; 


dll. 


co. 


+  (u  +  dx)(ox  4- 

C  X  cx 


which  reduces  to, 

dU'cx  “  “  [  77  {uox)\dxdydz. 


(25) 


Continuing  with  the  same  procedure  for  the  other  components  of  shear  stresses;  the  y- 
direction  and  the  z-direction,  the  total  change  in  work  due  to  normal  and  shear  stresses 
can  be  written  as  shown  in  equation  (26). 

dlV- 

I  “N  | 

-dt  1  —  (uax  +  IT^.  +  vvtjj,)  +  -j-  (uzyx  +  vcry  -I-  wxy2)  +  -j-  (ut2x  +  vr^.  +  wo2)  (26) 


The  total  energy  per  unit  mass  within  the  control  volume  is  the  sum  of  the  internal  and 
kinetic  energies,  given  by, 


Total  Energy  =  e  + 


(27) 


The  variation  in  kinetic  and  internal  energy  for  the  control  volume  is  shown  in  equations 
(28)  and  (2Vj,  respectively. 

d^iniernal  ~  d{pe)dxdydz  (2S) 
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(29) 


V2 

dE kinetic ^  d(p  2  )dxdydz 

Rewriting  these  two  equations  and  summing  them,  the  variation  of  total  energy  can  be 
written  as; 


dE 

di 


D/ 


V 

(pe  +  p  -t—  )dxdydz. 


(30) 


Changes  in  heat  content  due  to  conduction  only  are  considered.  According  to  Fourier’s 
Law  the  heat  flux  is  proportional  to  the  temperature  gradient,  so  that  the  heat  variation 
due  to  conduction  can  be  written  as; 


J_dQ__ 

A  di  q 


(31) 


Thus,  by  equating  the  amount  of  heat  transferred  into  the  volume  with  the  amount  of 
heat  leaving  the  volume,  the  following  relation  is  obtained; 

—  k  dvd:  +  (k  -~r~  k  -rr~  dx)dydz 
cx  '  cx  cx 

which  gives  the  heat  flux  in  the  x-direction, 

P  cT 

k  dxdvdz.  (32) 

Z'  V  z~-  V*  *  '  ' 


Repeating  similar  procedures  for  the  y-direction  and  the  z-direction  a  final  expression  for 
the  total  heat  variation  is  given  by; 


dQ_ 

dl 


cT  s  dT  \  c  /, 

-7— )  +  c(k  -7—  )  +  -7—  (k 

CX  CT  CZ 


)  dxdvdz. 
CZ  '  J 


(33) 


Substituting  equations  (26),  (30)  and  (33)  back  into  equation  (24),  the  general  energy 
equation  is  obtained.  [Ref.  3,19] 


Di  2  or  v  c x  '  cy  cy  '  cz  K  cz  ’ 


+  77  (uox  +  VTxy  +  »T„)  +  -j-  (UTyx  +  va  +  ttv7)  +  —  (UT„  +  VT„  +  wo,)  (34) 

C  Cl  CZ 
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B.  CONSERVATION  LAW  FORMULATION 

The  primitive  variable  formulation  of  the  Navier-Stokes  equations  shown  in  the 
previous  section  can  be  put  into  conservation  law  form  using  vector  identities.  The 
conservation  law  form  can  be  also  derived  by  applying  conservation  principles  on  a 
control  volume.  Because  physical  insight  is  gained  by  this  procedure  the  conservation 
law  form  derivation  is  outlined  in  the  next  section. 

This  formulation  stems  from  the  fact  that  certain  quantities  (i.e.  mass,  momentum, 
energy)  for  a  fluid  in  motion  are  conserved.  Conservation  implies  that  the  flux  of  a 
quantity  crossing  a  control  surface  and  the  net  effect  of  internal  sources  results  in  a 
variation  of  the  conserved  quantity.  These  sources  and  fluxes  depend  on  time  and  space 
as  well  as  fluid  motion.  The  fluxes  are  vectors  for  a  scalar  quantity  and  tensors  for  a 
vector  quantity.  Mass  and  energy  are  examples  of  a  scalar  quantity  whereas  momentum 
is  a  vector  quantity.  Molecular  motion  and  convective  transport  of  a  fluid  contribute 
to  flux.  Molecular  motion  has  the  tendency  to  make  the  fluid  homogeneous  and  has  a 
diffusive  effect. 

1.  General  Form  of  Conservation  Law 
a.  Scalar  Conservation  Law 

Considering  a  scalar  quantity  V  within  a  control  volume  V  ,  the  time  vari¬ 
ation  of  the  quantity  i'  is; 


This  should  be  equal  to  the  incoming  fluxes  (F  =  U 10  through  a  surface  S  (where  n  is 
the  unit  normal  vector  pointing  outward), 

-  j5«  •  FdS  =  —  J5F  .  is 

plus  any  possible  contribution  from  sources  of  U  . 
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Figure  5.  Flux  Diagram 

The  flux  vector  F  has  two  components,  a  diffusive  contribution  and  a  convective  con¬ 
tribution.  The  sources  can  be  written  as  the  addition  of  volume  sources  Qv  and  surface 
sources  Qs  ; 


\vQ>M’  +  \sQs-dS 

so  that  the  final  conservation  equation  for  the  scalar  quantity  U  is, 

jj ),  Wl'  =  ]fiydV  +  f5Cs  •  dS  -  J/  .  is.  (35) 

When  using  Gauss’s  Theorem,  equation  (35)  can  be  rewritten  as; 

I  - J„( IvtV  +  J,V . QsdV. 
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For  an  arbitrary-  volume  V  the  differential  form  of  the  conservation  law  is  given  by 
equation  (36). 


-jj-  +  V-(F~Qs)  =  Qv  (36) 

b.  Vector  Conservation  Law 

As  stated  earlier  if  the  conserved  quantity  U  is  a  vector,  then  the  flux  and 
surface  source  become  tensors;  F  =  U  •  V ,  Fs  ,  and  the  volume  source  in  turn  becomes 
a  vector,  Qv  .  An  analogous  derivation  as  in  the  conservation  of  a  scalar  quantity  can 
be  done  for  a  vector  quantity,  whose  integral  and  differential  form  are  shown  below. 

-J-  \tdV  +  jsF  .  £  =  jyQydV  +  !sQs  •  dS  (37) 

■%l  +  ^-(F-Qs)  =  Qv  (38) 

In  equation  (37)  the  convective  component  of  the  flux  tensor  can  be  written  in  tensor 
form  as; 


where  v  is  the  velocity  vector.  The  diffusive  component  of  the  flux  for  a  homogeneous 
system  can  be  written  as; 


where  k  is  the  diffusivitv  constant.  Equation  (35)  or  (36)  is  the  basic  formulation  of  the 
conservation  law  for  a  general  case.  When  continuity  of  flow  properties  is  assumed  (i.e. 
no  shocks  present),  then  equations  (36)  and  (38)  are  valid.  [Ref.  5:  pp.  25-55] 

2.  Equation  of  Mass  Conservation 

In  this  particular  instance  the  property  U  is  mass  and  no  diffusive  flux  is  pres¬ 
ent,  only  convection.  Therefore  equation  (35)  can  be  directly  written  as; 

or  in  differential  form  as  in  equation  (39).  [Ref.  5:  p.  33] 
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(39) 


-fr  +  V-(pV)  =  0 

3.  Equation  of  Momentum  Conservation 
For  this  case  the  conserved  quantity  is  momentum  which  is  a  vector.  From 

Newton's  Second  Law  it  was  mentioned  that  change  of  momentum  is  due  to  external 
volume  forces  and  internal  forces.  Assuming  a  Newtonian  fluid,  the  stresses  can  be 
written  as; 

a  =  —  pi  +  r 

a* 

where  1  is  the  unit  tensor,  so  that  —pi  is  the  hydrodynamic  pressure  along  the  diagonal. 
The  t  term  is  the  viscous  shear  stress  tensor,  equation  (15),  which  is  written  as; 

T/7  =  M((^;  +  5//)-y(V.n^). 

Referring  to  equation  (37)  and  assuming  that  the  external  volume  forces  is  zero  the  in¬ 
tegral  form  for  the  conservation  of  momentum  is; 

and  applying  Gauss's  Theorem, 

Jr4-  PVdV  +  \V  .  (p  V®  V)dV  =  .  adV 

where  ®  indicates  the  tensor  product  of  two  vectors.  This  can  be  written  in  differential 
form  as  shown  below.  [Ref.  5:  pp.  40-50.] 

4-{pV)  +  V»(pV®i'  +  pI  -t)«0  (40) 

4,  Equation  of  Energy  Conservation 

The  quantity  being  conserved  is  energy,  £  and  from  the  First  Law  of 
Thermodynamics  the  variation  in  energy  must  balance  with  the  work  of  the  forces  acting 
on  the  system  including  any  heat  addition.  The  convective  flux  of  energy  can  then  be 
written  as; 

Fc  =  pVE 


IS 


where  £  is  the  sum  of  the  internal  energy  plus  kinetic  energy.  Definition  of  the  diffusive 
flux  term  describes  that  diffusion  of  heat  for  a  fluid  at  rest  is  due  to  molecular  thermal 
conduction,  and  using  Fourier's  Law  of  Heat  Conduction,  the  diffusive  flux  term  can 
be  written  as; 

FD  =  kV.VT 

where  k  is  the  thermal  conductivity  and  T  is  the  absolute  temperature.  Assuming  no 
radiation,  chemical  reactions  or  work  due  to  external  forces,  again  the  Qv  volume  source 
is  zero.  The  net  work  done  on  the  fluid  by  the  internal  shear  stresses  acting  on  the  sur¬ 
face  of  the  control  volume  is  given  by; 

Qs-o.V. 

Using  equation  (37)  and  substituting  the  quantities  obtained  above,  the  equation  for 
energy  conservation  can  be  written  as; 

4- 1  \fiEdV  +  \spEV  .dS  =  \skTV.dS  +  \sCo.V).  dS 

or  in  differential  form  as  in  equation  (41). 

4-{pE)  +  V  *{V{Ep  +  p)-kVT- t  .  K)-0  (41) 


Equation  (41)  can  be  rewritten  as; 


P^~  +  PV  .V  =  4-(k-~)  +  4-(k-^)  +  4-(k4r~)  +  pO 
Di  cx  K  cx  '  cy  cy  cz  cz' 


where  O  is  called  the  dissipation  function.  Dissipation  represents  the  heat  equivalent 
of  the  rate  at  which  the  mechanical  energy  is  lost  during  deformation  of  the  medium  due 
to  viscosity.  The  dissipation  function  is  given  by; 


o  =  2r(4^)2  +  (4-)2  +  (4^)2]  +  (4-+-^-)2  + 

_cx  cy  c:  J  cx  cy 


)>  +  ( i!L  +  iH- )>  _  2  ( 4"  +  it  + 

c:  cx  3  cx  cv  cz 
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The  system  of  partial  differential  equations  given  by  equations  (39),  (40)  and  (41)  can 
be  written  in  compact  vector  notation  as; 

dq  —  — 

-T-  +  V  .  Q  =  0  (42) 

cl 

where  q  is  the  vector  of  dependent  conservative  variables  and  Q  is  a  vector  composed 
of  the  nonlinear  inviscid  and  viscous  fluxes.[Ref.  5:  pp.  45-50] 

5.  Strong  Conservation  Form 

The  strong  conservation  law  form  given  by  equations  (39),  (40)  and  (41)  in 
vector  notation  can  be  written  for  a  Cartesian  coordinate  system  as; 


eg  cl: '  cF  cG  BR.dS  ST 

-i  I  1  <■*  " '  I  I  ^ 

ct  c.x  cy  oz  cx  cy  cz 


(43) 


where  q  is  the  vector  of  conservative  variables  and  E,  F  ,  and  G  are  the  flux  vectors 
given  by; 


g  = 


P 

pu 

pv 

pw 

e 


(44) 


pu 

pv 

pw 

pu2  +p 

pvu 

pwu 

puv 

F  = 

pv2  +  p 

G  = 

pwv 

pun- 

pvw 

pw2  +  p 

ip  +  e)u 

ip  +  e)v 

{p  +  e)w 

(45) 


The  vectors  of  R,  S  ,  and  T  contain  the  viscous  terms.  When  they  are  omitted,  the  Euler 
equations  are  recovered. 
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(46) 


0 

0 

0 

rxx 

T>* 

r*y 

S  = 

Tyy 

T  = 

x*y 

*X2 

_  V 

^22 

(f-P  +  qc)^ 

1 

+ 

1^ 

• 

1  h-* 

(T  •  V  +  qc)2_ 

The  product  term  of  T  •  V  is  written  in  component  form  as  below. 

( T  •  V)x  =  ?xxu  +  ?xyv  +  tx:w 

(  T  •  V)y  =  T yXU  +  Tyy  v  +  t>2w  (47) 

o 

(T'  V)2  =  x2Xu  +  x2yv  +  x22w 

The  heat  flux  vector  qc  is  the  heat  transfer  by  conduction  and  can  be  written  as; 

qc  =  -kVT=-K(ala^a22)T  (48) 

where, 


K  = 


n 

P>{y  -  1) 


In  the  above  equations  a  denotes  the  speed  of  sound,  Pr  is  the  Prandtl  number,  cp  is  the 
specific  heat  at  a  constant  pressure  and  e  is  the  total  energy  per  unit  volume.  [Ref.  20] 
Pressure  and  energy  are  related  by  the  perfect  gas  law  as  follows. 

p  ~(y  -  1)£*  -■y  («z  +  v2  +  w2)j 

These  equations  can  be  transformed  into  different  curvilinear  coordinate  systems  in  or¬ 
der  to  facilitate  the  numerical  implementation. 

A  coordinate  mapping  is  introduced  which  allows  the  transformation  of  the 
equations  of  motion  from  a  Cartesian  coordinate,  time  varying,  nonorthogonal  coordi¬ 
nate  system.  The  mapping  is  linked  to  the  Cartesian  coordinates  as  follows; 

(  =  <f  (XJ,2.l) 
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n  =  n{xy,2,t) 


C  =  «w./). 


The  Cartesian  coordinate  system  is  the  physical  domain,  and  the  transformed  space  is 
referred  to  as  the  computational  domain.  This  computational  domain  is  orthogonal  with 
a  uniform  rectangular  mesh  so  that  unweighted  differences  can  be  taken  to  form  the 
derivatives. 

The  thin  layer  compressible  Navier-Stokes  equations  are  obtained  from 
equation  (43)  by  retaining  the  viscous  terms  only  along  the  direction  that  is  normal  to 
the  body.  Also,  the  derivatives  of  the  stress  terms  in  the  crossflow  (i.e.  y,  z )  directions 
are  discarded.  The  thin  layer  formulation  of  the  strong  conservation  law  form  of  the 
governing  equations  for  a  curvilinear  coordinate  system  (£,  tj,  £)  along  the  axial, 
circumferential,  and  normal  direction,  respectively  can  be  written  as; 


where  q .  F,  G,  H  ,  and  S  are, 


<7 


AAA  A 


cF  cG  dli  _  1  cS 

T  v  i  «■»  i  -n  v  n  v 


cc 

*+-  ~  "T  “ 

C)J 

Re  e; 

P 

pU 

pu 

pul'  +  ZxP 

pv 

A  1 
fa7 

PVU  +  tyP 

pw 

pwl'  + 

e 

(e  +  p)U-Z,p 

(49) 


pV 

pJV 

puV+tij 

pu\V  +  Ijj? 

p\V+VyP 

prfV+ZyP 

pwV+iij> 

pwW+lj) 

{e  +  p)V-q,p 

-  — 

{e  +  pW-tf 
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0 


s 


Atmjt/j  +  (M/3)m2Cx 

nmxv{  +  (#i/3)/n2Cy 

MW,w{  +  (m/3)w2£, 
nmxm3  +  (Ai/3)m  +  2(Cxw  +  +  C7*') 


Furthermore,  it  is  defined  that; 

"h  =  c*  +  c2  +  c2 

"»2  =  +  C7W{ 

A  2 

,  2  ,  2  ,  2Wa  ,  K  ,  Cfl  v 

m3  =  (u  +  v  +  w  )l2  +  y^(  -jr ) 

and  U,  V  ,  and  IV  are  the  contravaTiant  velocity  components  given  by, 

L  =  UC x  4-  \'ty  +  wt2  "b  *t  r 

f/=u^  +  ’^  +  Mlz  +  tlt 

w=u:x  +  v:y  +  w:2  +  c,. 

Again  analogous  to  the  previous  derivations  the  pressure  is  related  to  density  and  total 
energy  through  the  equation  of  state  for  an  ideal  gas. 

C.  NUMERICAL  IMPLEMENTATION 
l.  The  Numerical  Algorithm 

The  solutions  over  a  strake-delta  wing  configuration  resembling  a  modern 
fighter  aircraft  planform  will  be  presented  in  the  last  part  of  this  thesis.  Even  though  the 
main  effort  of  this  work  was  not  the  numerical  solution  of  the  governing  equations  (i.e. 
the  compressible  Navier-Stokes  equations),  the  technique  used  for  the  numerical  imple¬ 
mentation  is  briefly  described  in  the  following  paragraphs. 

The  numerical  scheme  used  for  the  solution  of  the  governing  equations  is  based 
on  a  finite  difference  discretization  of  the  thin  layer  Navier-Stokes  equation  [Ref.  6). 
The  numerical  integration  was  performed  using  a  partially  flux-split  numerical  scheme. 
L’pwinding  was  performed  in  the  main  flow  direction  using  flux  vector  splitting  while 
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central  differencing  was  used  in  the  other  two  directions.  The  factored  form  of  the  re¬ 
sulting  algorithm  is  as  follows; 


[  /  +  hb\{A+)n  -l-  h6iCn  -  hRe~]6:r]MnJ  -  Z>,|{] 
x  [/  +  hlf>{A-)n  +  h6nB*  -  Dt  = 

-  Atff  [(F*)"  -  C]  +  ^nn  ~  Kol  +  W  -  GJ  +  6 £Hn  +  b£Sn  -  SJ) 


~  De(qn  -  O- 


The  explicit  dissipation  D,  was  used  along  the  directions  where  central  differencing  was 
applied.  The  implicit  dissipation  term  Z),  was  added  for  numerical  stability.  Steady 
aerodynamic  flows  at  subsonic  flows  (M  =  0.2)  do  not  contain  shock  waves  and  can 
be  quite  well  predicted  by  a  central  difference  scheme  that  is  augmented  by  these  dissi¬ 
pation  terms. [Ref.  20  ] 

2.  Turbulence  Model 

Simulation  of  high  Reynolds  number  flows  is  obtained  by  the  solution  of  the 
Reynolds  averaged  Navier-Stokes  equations.  These  equations  have  extra  unknowns  and 
are  commonly  called  the  Reynolds  stresses  [Ref.  3).  The  relations  between  the  Reynolds 
stresses  and  the  mean  flow  quantities  is  the  well  known  closure  problem.  In  practice 
some  turbulence  model  is  used  which  relates  the  Reynolds  stresses  with  the  mean  flow 
quantities.  The  turbulence  model  selected  for  this  research  was  an  algebraic  eddy 
viscosity.  This  model  is  the  Baldwin- Lomax  model  as  modified  by  Degani  and  Schiff  to 
treat  three-dimensional  separated  flows  [Ref.  21,22). 

The  turbulence  is  simulated  in  terms  of  an  eddy  viscosity  coefficient  n,  .  The 
coefficient  of  viscosity  and  the  heat  flux  term  in  the  Navier-Stokes  equations  are  re- 

n  n 

placed  with  m  +  M,  and  —  +  —  ,  respectively.  The  turbulence  model  is  similai  to  one 
developed  by  Cebeci  with  modifications  that  allow  for  the  locating  of  the  boundary  layer 
[Ref  23).  A  two  layer  algebraic  eddy  viscosity  model  is  used  where  the  Prandtl-Van 
Driest  formulation  is  used  in  the  inner  region  and  the  Clauser  formulation  is  used  in  the 
outer  region  [Ref.  21).  The  inner  region  is  any  normal  distance  from  the  wall,  y  ,  that 
is  less  than  or  equal  toy,™....- .  If  this  is  the  case  then  n,  is  defined  by  the  following  ex¬ 
pression; 

(burner  =  I  «*>  I 
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where, 


l=ky 


1  -  exp(  — ) 
A 


and, 


i  ./  cu  cv  ^  ,  ,  cv  cm-  a  ,  ,  cw  cu  ,2 

w  =  /( ~ - —  )  +(— - T—  )  4-  ( — - 7-  ) 

v'  cy  c  x  cz  cy  cx  cz 


+ 

y 


p*u-y 


Jp*r*y 

Pw 


Ify  is  greater  than yer,„m,r  then  n,  is  defined  by; 

(Pt)  outer  ~  KCcpF^/a^gF Klebi}) 


where  K  is  the  Clauser  constant,  Cc,  is  an  additional  constant,  and  for  boundary  layers, 

F\\'ake  ~  Tmax^max 


or  for  wakes  and  separated  boundary  layers, 


Fll'akc  ~  C frji) n 


1 

L  Dif 
f 

1  max 


In  the  above  equation  UDl,  is  the  difference  between  the  maximum  velocity  atym„  and 
the  minimum  velocity  in  the  profile.  The  quantities  of ymtx  and  Fm,x  are  calculated  using; 


F(y)  =y  |  w  1 


_  .+ 

1  -  exp(  — ~ ) 
A 


The  function  FKlti(y)  is  the  KlebanofTintermittency  factor  and  is  defined  as; 

All  other  parameters  are  constants  determined  empirically  and  given  in  [Ref  21J.  The 
use  of  this  model  eliminates  the  need  for  finding  the  edge  of  the  boundary  layer  and  re¬ 
duces  one  of  the  sources  of  error  in  the  Navier-Stokes  solutions. 
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III.  SURFACE  GRID  GENERATION 


The  focal  point  of  this  work  was  the  generation  of  the  computational  mesh  over  a 
three-dimensional  strake-delta  wing  configuration  that  models  a  modem  fighter  aircraft 
planform.  Therefore,  in  the  next  two  chapters  the  surface  and  field  grid  generation 
procedures  are  described.  The  considerations  which  must  be  taken  into  account  in  order 
to  construct  a  surface  grid  suitable  for  the  subsequent  generation  of  the  field  are  dis¬ 
cussed.  Finally,  the  various  field  grid  generation  methods  are  discussed  and  two  different 
approaches  which  were  used  to  generate  the  field  grid  over  the  strake-delta  wing  con¬ 
figuration  are  described. 

The  first  step  to  be  taken  in  establishing  a  finite  difference  or  finite  element  scheme 
for  solving  a  system  of  partial  differential  equations  is  to  replace  the  continuous  domain 
by  a  finite  mesh,  commonly  known  as  a  grid.  Grid  generation  is  one  of  the  central 
problems  in  the  procedure  to  obtain  a  numerical  solution.  A  well  constructed  grid 
greatly  facilitates  the  numerical  solution  of  a  system  of  P.D.E.s.  On  the  other  hand,  an 
improper  grid  choice  may  lead  to  instabilities,  inaccuracies  and  or  lack  of  convergence. 
Numerical  grid  generation  is  a  procedure  for  the  orderly  distribution  of  observers  over 
the  physical  field  domain  in  such  a  way  that  efficient  communication  among  the  ob¬ 
servers  is  possible.  Also,  it  assures  that  all  physical  phenomena  of  interest  in  the  entire 
field  may  be  represented  with  sufficient  accuracy  by  this  finite  collection  of  observers. 

Grid  generation  for  two-dimensional  domains  is  relatively  simple  and  may  be 
achieved  with  purely  algebraic  techniques,  even  for  relatively  complex  domains  [Ref. 
24).  In  addition,  for  suitable  geometries  of  the  boundaries  conformal  mapping  techniques 
may  be  used.  Conformal  mapping  techniques  have  the  advantage  that  they  are  relatively 
simple  and  inexpensive  [Ref.  25  :  pp.  48S-490][Ref.  4:  pp.  7-56].  They  also  preserve  grid 
orthogonality,  but  their  use  is  limited  to  domains  with  simple  boundaries  where  a  con¬ 
formal  transformation  between  the  physical  domain  and  a  simpler  transformed  domain 
may  be  readily  defined. 

The  generation  of  a  computational  grid  for  a  three-dimensional  domain,  however, 
presents  greater  difficulties.  For  a  limited  class  of  external  and  internal  domains  it  is 
sometimes  possible  to  fill  the  entire  three  dimensional  domain  with  a  sequence  of  two- 
dimensional  plane  grids  that  will  constitute  the  entire  three-dimensional  field  grid.  An 
application  of  this  idea  is  shown  later  for  the  construction  of  the  field  grid  over  the 
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double-delta  wing.  In  many  instances  however,  it  is  either  difficult  to  decompose  the 
three-dimensional  domain  into  a  sequence  of  two-dimensional  domains,  or  it  is  prefera¬ 
ble  to  construct  a  purely  three-dimensional  mesh.  Of  course,  both  the  complexity  and 
computing  time  of  a  three-dimensional  grid  generation  method  will  be  higher.  In  any 
case,  the  definition  of  the  surface  boundaries  must  be  done  precisely  and  accurately. 

Before  any  work  can  be  started  on  a  three-dimensional  field  grid  over  a  body  one 
must  first  define  the  surface  geometry  of  the  body.  The  definition  of  the  body's  surface 
and  its  quality  is  imperative  to  the  success  of  the  field  grid.  There  exist  many  avenues 
to  create  the  surface  grid,  some  include  algebraic  techniques,  cubic  Hermite  functions, 
Bezier  curves  or  Non-Uniform  Rational  B-splines  (NURBS)  [Ref.  4:  pp.  237-249]  [Ref. 
25:  pp.  497-503].  Each  of  these  techniques  has  its  advantages  and  the  exact  method  that 
will  best  fit  a  particular  surface  will  vary.  The  availability  of  accurate  data  for  de¬ 
scription  of  the  surface  geometry  will  also  play  a  major  role  in  the  generation  of  the 
surface  grid.  If  all  surfaces  can  be  defined  in  terms  of  equations,  then  an  algebraic 
technique  might  prove  to  be  the  most  efficient.  Whereas,  if  the  surface  is  very  complex, 
as  is  the  case  for  actual  aircraft  surfaces,  all  that  is  available  are  two-dimensional  cross- 
sections  and  a  curve  fitting  technique  will  have  to  be  used.  Whatever  the  method,  it  is 
of  utmost  importance  that  the  surface  grid  generation  program  be  written  in  such  a  way 
that  it  will  enable  maximum  flexibility  in  the  number  of  grid  points  and  their  distrib¬ 
ution.  This  early  concern  and  respect  for  versatility  will  pay  large  dividends  upon  sub¬ 
sequent  generation  of  the  field  grid.  For  even  after  the  surface  grid  has  been  completed, 
an  interactive  trial  and  error  process  of  changing  the  surface  grid  will  be  required  to 
achieve  an  effective  field  grid. 

A.  DOUBLE-DELTA  WING  SURFACE  GRID 

The  dimensions  of  the  double-delta  wing  model  are  shown  in  Figure  (7).  From  this 
figure  it  can  be  easily  seen  that  most  of  the  surfaces  can  be  defined  by  linear  relation¬ 
ships  with  the  only  exception  being  the  NACA  64-005  cross-section.  For  this  reason 
algebraic  grid  generation  on  the  surface  was  chosen.  For  a  linear  relationship  and  ana¬ 
lytically  defined  points  no  advantage  is  gained  by  using  a  curve  fitting  method.  The  part 
of  the  wing  that  contains  the  NACA  64-005  airfoil  cross-section  required  special  treat¬ 
ment.  Generation  of  the  surface  grid  over  the  airfoil  cross-section  would  require  a  curve 
fitting  technique.  Much  work  has  been  done  in  this  area  and  NURBS  can  produce  ex¬ 
cellent  results  in  approximating  airfoils.  The  main  advantage  of  this  technique  is  that 
it  provides  the  flexibility  of  modifying  the  cross-section  or  shape  of  the  surface  by  simple 
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changes  of  user  specified  parameters.  The  disadvantage  is  that  the  complexity  is  higher 
and  the  redistribution  of  the  points  that  represent  the  airfoil  is  also  more  difficult.  Be¬ 
cause  the  cross-section  is  a  NACA  airfoil  whose  contour  shape  can  be  well  approxi¬ 
mated  by  straight  lines,  the  use  of  a  purely  algebraic  technique  for  the  entire 
double-delta  wing  surface  grid  was  utilized. 

Once  this  decision  was  made,  the  areas  containing  singularities  had  to  be  identified 
and  the  program  for  the  algebraic  grid  generation  had  to  be  written.  Because  the  surface 
grid  is  used  as  initial  or  boundary  conditions  for  the  generation  of  the  field  grid,  much 
care  had  to  be  taken  to  avoid  any  singularities  that  would  propagate  into  the  field  grid. 
In  addition,  special  care  must  be  taken  at  the  regions  of  sharp  comers  such  as  the  lead¬ 
ing  edge.  In  these  areas  it  is  not  possible  to  maintain  field  grid  orthogonality  and  the 
location  of  these  acute  angles  can  be  seen  in  Figure  (7).  In  general,  these  areas  are  found 
on  the  entire  leading  edge,  at  the  apex  and  at  the  rectangular  edge  near  the  wingtip. 
These  corners  had  to  be  approximated  by  "rounding"  ofT  these  areas  with  a  radius  that 
was  very  small.  The  radius  used  was  0.001%  of  the  chord  length,  so  to  the  naked  eye 
the  surface  grid  appears  to  be  a  sharp  corner.  This  rounding  of  acute  angles  allows  the 
field  grid  to  maintain  orthogonality  which  is  a  desirable  feature  for  subsequent  numerical 
implementation;  see  Figures  (11),  (14),  (17),  (20)  and  (22).  Also,  high  grid  resolution  is 
provided  at  the  same  time  in  these  areas  where  the  change  of  the  flow  field  variables  is 
expected  to  be  rapid.  The  methodology  for  generating  the  source  code  that  would 
compute  the  grid  points  was  to  progress  from  the  nose  in  an  axial  direction  through  to 
the  wake.  The  grid  points  were  essentially  generated  for  a  two-dimensional  cross-section 
in  the  yz-plane,  then  an  incremental  step  in  the  x-direction  was  made  and  again  the  grid 
points  for  a  new  yz-plane  were  computed.  The  source  code  was  written  in  five  logical 
sections  that  defined  regions  of  the  wing  with  similar  cross-sections.  These  five  sections 
were  the  apex,  the  strake,  the  wing,  the  trailing  edge  rectangular  section  and  the  wake. 

1.  The  Apex 

Special  care  had  to  be  taken  in  the  modeling  of  the  nose  region.  The  apex  of 
the  wing  is  a  single  point  that  transitions  to  a  diamond  shape  cross-section.  Taking  into 
account  that  smoothness  has  to  be  maintained,  a  hemisphere  may  be  used  to  provide 
smooth  transition  between  the  singular  point  of  the  apex  and  the  diamond  cross-sections 
at  the  nose,  see  Figure  (11).  The  radius  of  this  hemisphere  is  0.001  %  of  the  chord  length 
and  allows  for  a  smooth  transition.  The  radius  of  the  sphere,  the  number  of  grid  points 
for  the  axial  (x-direction)  and  the  circumferential  (y-direction)  were  inputted  by  the  user. 
An  incremental  angle  was  determined  for  both  the  xz-planes  and  yz-planes  which  then 
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enabled  the  computation  of  the  grid  points.  The  yz-plane  cross-sectional  grid  points 
were  then  calculated  using  incremental  yz-plane  angles  as  the  x-location  progressed 
downstream  using  the  incremental  xz-plane  angles. 

2.  The  Strake 

The  main  concern  in  defming  this  section  was  that  the  leading  edge  has  a  corner 
that  is  relatively  sharp  and  would  preclude  a  field  grid  that  is  orthogonal.  Therefore,  the 
sharp  leading  edge  was  approximately  rounded  as  shown  in  Figure  (12)  through  Figure 
(17).  The  computations  here  involved  a  user  specified  radius  that  was  the  same  as  the 
one  in  the  approximation  of  the  apex.  This  radius  was  maintained  to  allow  a  smooth 
transition  from  the  sphere  of  the  apex  to  the  diamond  cross-section  of  the  strake.  The 
surface  grid  generator  code  provides  the  versatility  to  change  the  number  of  grid  lines  in 
this  radius  which  is  kept  constant  for  every  cross-section.  This  issue  becomes  important 
as  the  ratio  of  radius  to  wing  span  drastically  changes  between  the  apex  and  junction 
of  the  strake  and  wing.  Again  similar  logic  to  the  one  used  to  define  the  apex  was  used 
here.  The  difference  being  that  the  increment  of  the  x-coordinate  was  computed  de¬ 
pending  on  the  number  of  grid  lines  along  the  x-direction  used  to  define  the  strake  part 
of  the  body.  Simple  relations  from  analytic  geometry  are  used  to  determine  the  y  and 
z-coordinate  as  a  function  of  the  x-location. 

The  distribution  or  clustering  of  the  grid  lines  will  be  discussed  in  more  detail 
later  in  this  section.  It  is  important  to  mention  that  the  distances  between  successive 
grid  points,  along  the  x  and  y-directicns  was  determined  by  calling  a  subroutine.  This 
allows  to  experiment  with  many  different  distributions  with  a  simple  change  of  input 
parameters. 

3.  The  Wing 

For  the  wing  section  three  areas  required  special  attention.  The  first  was  like 
the  strake.  in  that  the  leading  edge  forms  a  comer  which  unlike  the  strake  was  not  as 
sharp.  This  was  because  of  the  NACA  64-005  cross-section  has  a  fmite  curvature  at  the 
leading  edge.  Part  of  the  leading  edge  did  not  require  rounding  and  the  number  of  grid 
lines  approximating  the  edge  was  reduced,  see  Figure  (15)  through  Figure  (17).  The 
NACA  64-005  cross-section  is  generated  by  a  subroutine  which  requires  as  input  only 
the  normalized  root  chord  length  of  the  airfoil,  xe  =  jct(y)  •  The  area  of  the  wing  spanning 
between  the  wing  centerline  and  the  part  of  the  wing  having  a  NACA  64-005  cross- 
section  was  defined  by  linear  interpolation.  Some  sort  of  curve  fitting  method  could 
have  been  used,  but  since  the  wing  is. thin  and  the  distance  is  short,  a  linear  approxi¬ 
mation  was  assumed  to  be  sufficiently  accurate. 
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4.  The  Trailing  Edge  Rectangular  Section 

This  particular  section  required  the  most  challenging  surface  grid  definition  in 
the  early  stages.  This  was  primarily  because  the  wingtip  has  a  variable  finite  thickness 
that  depends  on  the  x -location,  see  Figure  (18)  through  Figure  (20).  The  actual  calcu¬ 
lation  of  the  grid  point  locations  was  simple,  but  the  number  of  grid  lines  at  the  wingtip 
region  had  to  change  due  to  this  variable  thickness.  This  required  to  take  grid  lines  out 
of  the  edge  and  redistribute  them  back  onto  the  upper  and  lower  surfaces.  This  is  the 
reason  why  the  grid  lines  in  this  section,  when  viewed  in  the  xy-plane  appear  staggered, 
see  Figure  (9).  The  trailing  edge  has  a  finite  thickness  (0.002%  of  chord  length)  and  thus 
avoids  any  singularities  that  unnecessarily  complicate  the  subsequent  numerical  imple¬ 
mentation  and  the  generation  of  the  field  grid. 

5.  The  Wake 

For  the  numerical  implementation,  an  extension  of  the  far  end  of  the  computa¬ 
tional  domain  of  2.0  -  3.0  root  chords  beyond  the  body  is  required.  The  wake  was  rela¬ 
tively  simple  to  generate  because  all  that  varied  was  the  x-coordinate.  All  the  yz-planes 
remain  constant  from  the  trailing  edge  to  the  end  of  the  grid.  Examples  of  this  cross- 
section  can  be  seen  in  Figure  (21)  and  Figure  (22).  The  wake  extends  for  2.0  root  chord 
lengths  beyond  the  wing  trailing  edge.  This  length  was  selected  because  it  allowed  for 
a  smooth  transition  from  the  wing  to  the  wake  for  a  given  number  of  x-direction  grid 
lines. 

B.  DISTRIBUTION  PARAMETERS 

Flexibility  in  the  distribution  of  the  surface  grid  points  in  both  the  x-direction  and 
the  y-direction,  which  are  shown  in  Figure  (10),  is  important  for  the  surface  grid.  The 
foremost  problem  is  to  ensure  that  the  distance  between  successive  grid  points  makes  a 
smooth  transition.  Of  course,  the  spacing  between  the  surface  grid  points  could  have 
been  made  the  same,  but  this  is  impractical  because  the  total  number  of  grid  lines  would 
be  excessive  due  to  the  small  radius  that  was  used  to  approximate  the  corners.  Therefore 
a  distribution  of  grid  points  must  be  developed  that  is  very  dense  at  the  comers  and 
sparser  in  the  other  regions.  High  grid  clustering  is  also  required  in  areas  where  steep 
gradients  in  the  flow-field  are  expected,  such  as  the  leading  edge  where  the  leading  edge 
vortices  appear.  The  distribution  in  the  y-direction  can  be  seen  in  Figure  (9)  and  Figure 
(10)  for  various  cross-sections  of  the  wing  and  the  distribution  for  the  x-direction  can 
be  seen  in  Figure  ( 10).  The  y-direction  has  a  high  grid  clustering  around  the  leading  edge 
which  becomes  sparser  near  the  centerline  of  the  wing.  This  was  the  general  procedure 
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followed  for  the  y*direction  distributions  in  all  cross-sections.  The  distribution  along  the 
x-direction  required  a  higher  density  in  the  nose  region  and  sparser  distribution  at  the 
area  near  the  end  of  the  wake.  A  subtle  change  to  a  higher  density  occurs  where  the 
wing  has  a  geometry  change  as  can  be  seen  in  Figure  (9). 

Stretching  of  the  grid  points  along  a  coordinate  direction  can  be  obtained  by  using 
simple  algebraic  functions  such  as  linear,  exponential  mapping  or  trigonometric  func¬ 
tions.  Use  of  these  functions  allows  for  a  smooth  transition  from  sparse  grid  densities 
to  high  grid  densities.  A  quadratic  function  was  first  attempted  but  this  resulted  in  a 
distribution  that  became  less  dense  too  fast  and  would  spread  the  grid  points  out  to  an 
excessive  amount  near  the  centerline.  Next  a  linear  stretching  function  was  used  and  this 
gave  much  better  results  but  did  not  allow  a  "smooth"  transition  from  the  high  grid 
density  region  to  the  region  with  sparser  grids.  This  effect  was  more  pronounced  along 
the  y-direction.  The  linear  function  allows  a  more  constant  distribution  over  the  whole 
wing  in  the  x-direction.  The  linear  equation  shown  below  was  used  for  the  linear 
stretching: 


where  the  user  specifies  the  parameter  c  depending  on  the  desired  degree  of  stretching. 
A  value  ofc  =  1.0  would  result  in  an  equidistance  spacing,  whereas  a  c  —  2.0  would  result 
in  a  high  clustering  of  grid  points  near  one  or  both  ends.  Difficulties  were  not  en¬ 
countered  in  the  x-direction  because  the  transition  from  the  low  to  high  density  distrib¬ 
utions  were  not  as  extreme  as  in  the  y-direction.  By  changing  the  linear  stretching 
parameter,  a  smooth  transition  from  more  to  less  dense  areas  was  achieved.  The  grid 
stretching  in  the  axial  direction  required  different  values  of  c  depending  on  the  wing 
section;  for  example,  c  =  1.005  was  used  prior  to  the  trailing  edge  and  c  =  1.205  after  the 
trailing  edge  of  the  wing.  The  effect  of  these  constants  can  be  seen  in  Figure  (9).  Each 
representative  cross-section  had  to  be  investigated  and  a  constant  assigned  that  allowed 
a  smooth  transition  from  one  section  to  another.  These  constants  were  determined  by 
a  trial  and  error  procedure,  but  experience  gained  by  many  iterations  expedited  the 
process.  From  the  many  iterations  for  the  surface  grid  alone,  an  appreciation  for  the 
flexibility  of  the  source  code  was  gained. 

To  resolve  the  y-direction  distribution,  the  stretching  function  first  attempted  was 
exponential  which  proved  to  be  inadequate.  The  exponential  stretching  is  obtained  by 
using  the  following  expression; 
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yk+ i  =  wl 

where  c  and  s  are  parameters  chosen  by  the  user  to  produce  a  desired  distribution.  The 
dimensions  of  the  radius  used  to  approximate  comers  are  so  small  compared  to  the 
characteristic  dimension  (i.e.  the  chord  length)  that  for  the  desired  number  of  grid  points 
in  the  y-direction  the  transition  did  not  occur  smoothly.  Finally,  a  sinusoidal  distrib¬ 
ution  produced  better  results.  The  equation  below  was  used  to  obtain  this  distribution; 

yk  =  c  sin  b 

here  c  is  a  user  specified  parameter  and  b  is  allowed  to  incrementally  change  over  a 
specified  range  of  angles  in  order  to  obtain  the  desired  section  of  the  sine  curve.  The 
distribution  produced  the  best  results  for  a  sine  curve  segment  from  0  to  45  degrees. 
Freedom  was  written  into  the  source  code  to  use  constants  to  finely  adjust  the  distrib¬ 
ution  but  were  not  requiied  because  of  sufficient  results  without  them. 

C.  PROGRAM  FEATURES  FOR  THE  SURFACE  GRID 

The  source  code  for  the  surface  grid  can  be  seen  in  Appendix  E.  The  main  concern 
in  the  construction  of  the  source  code  for  this  problem  was  to  give  the  author  maximum 
flexibility  in  the  generation  of  this  surface  grid.  Listed  below  are  some  of  the  features 
that  can  be  easily  changed  via  an  input  file. 

•  The  number  of  grid  points  in  the  x-direction  at  five  different  sections. 

•  The  number  of  total  grid  points  in  the  y-direction. 

•  The  radius  used  in  the  approximations  of  sharp  edges. 

•  The  number  of  grid  points  used  in  the  radius  for  the  corner  approximations. 

•  The  sweep  angles  of  the  strake  and  the  wing. 

•  Maximum  widths  of  the  strake  and  wing. 

•  Lengths  of  the  strake,  the  wing,  the  rectangular  section  and  the  wake. 

•  Distribution  constants  at  five  locations  in  the  x-direction  and  at  six  locations  in  the 
y-direction. 

This  flexibility  in  the  surface  grid  paid  a  major  dividend  in  the  subsequent  generation  of 
the  field  grid.  This  is  because  the  surface  and  field  grid  generation  is  an  interactive 
process  that  usually  requires  changes  in  the  surface  grid.  Another  feature  is  that  the 
distribution  functions  are  written  as  subroutines.  This  allowed  the  author  the  flexibility 
of  trying  different  functions  based  on  the  geometry'  and  desired  gradients  of  distribution. 
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Another  subroutine  is  the  calculation  of  the  grid  points  that  are  part  of  the  N'ACA 
64-005  airfoil  cross-section.  This  enabled  the  changing  of  this  cross-section  by  merely 
changing  two  lines  of  the  data  statement.  The  subroutines  that  are  included  in  Appendix 
E  are  linear,  quadratic,  exponential  and  sinusoidal.  Again  it  is  emphasized  that  the 
source  code  for  the  surface  grid  generation  must  be  as  versatile  as  possible. 

This  program  can  be  used  to  generate  a  surface  grid  for  a  wing  with  dimensions  that 
are  different  from  the  one  chosen  by  the  author.  But  because  great  care  needs  to  be 
taken  in  grid  point  distribution  a  change  in  the  dimensions  would  most  definitely  require 
an  adjustment  of  radius,  the  number  of  points,  distribution  constants  and  even  distrib¬ 
ution  functions.  A  close  examination  of  every  grid  constructed,  either  visually  or  com¬ 
putationally,  must  be  completed  to  ensure  that  no  intersection  of  the  grid  lines  occurs. 
The  program  that  was  originally  written  was  continually  revised  to  permit  an  adequate 
construction  of  the  field  grids.  Included  in  Appendix  F  is  the  final  program  that  was 
utilized  for  the  generation  of  the  surface  grid  for  the  final  spherical  field  grid  topology. 
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IV.  FIELD  GRID  GENERATION 


In  this  chapter  different  numerical  techniques  for  the  field  grid  generation  are  pre¬ 
sented.  Their  advantages  and  disadvantages  as  far  as  grid  generation  and  numerical 
implementation  are  explained.  Most  available  field  generation  techniques  require  as  in¬ 
put  a  surface  grid  which  must  be  constructed  before  the  field  grid  generation.  The  gen¬ 
eration  of  the  field  grid  is  directly  dependent  on  the  distribution  of  the  grid  points  on  the 
body  surface.  There  are  several  ways  to  generate  field  grids,  a  few  of  which  include 
methods  involving  solutions  of  Elliptic  Partial  Differential  Equations  (P.D.E.),  Parabolic 
P.D.E.s,  or  Hyperbolic  P.D.E.s.  For  a  limited  class  of  problems  algebraic  methods  can 
also  be  used.  In  the  following  paragraphs  the  various  grid  generation  methods  based 
on  the  solutions  of  P.D.E.s  will  be  described.  The  hyperbolic  method,  which  was  used 
to  generate  the  field  grid  over  the  double-delta  wing  will  be  explained  in  detail,  whereas 
only  a  brief  description  of  elliptic  and  parabolic  techniques  will  be  given. 

The  classification  of  P.D.E.s  into  hyperbolic,  parabolic  or  elliptic  type  is  obtained 
from  the  general  form  of  the  quasi-linear  second  order  P.D.E.,  given  by; 

*uxx  +  buxy  +  cu>v  +  dux  +  euy  +fu  =  g  ( 1 ) 

here  w  =  is  the  dependent  variable  and  the  coefficients  a,b,c,d,ef  and  g  are  func¬ 
tions  of  x  andy  .  The  type  of  equation  (1)  is  determined  from  the  sign  of  the  quantity 
b2  —  4ac  as  follows. 


b2  —  4 ac  <  0 

( Hyperbolic ) 

(2) 

b 2  —  4ac  =  0 

(Parabolic) 

(3) 

b2  —  4 ac  >  0 

(Elliptic) 

(4) 

Each  type  of  equation,  hyperbolic,  parabolic  or  elliptic  has  certain  characteristic  prop¬ 
erties  which  can  be  successfully  utilized  for  the  grid  generation  in  two  and  three- 
dimensional  domains  [Ref.  4:  pp.  188-277].  For  example,  the  solution  of  an  elliptic 
equation  in  the  interior  of  a  domain  depends  on  the  specification  of  data  o\cr  the  entire 
boundary.  Therefore,  when  a  grid  is  generated  by  the  solution  of  an  elliptic  P.D.E.  all 
the  boundary  data  must  be  specified. 
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The  main  feature  of  a  hyperbolic  type  of  problem  is  that  the  solution  starts  from  an 
initial  condition  and  propagates  in  time  along  certain  directions  known  as  the  charac¬ 
teristic  directions.  Utilization  of  this  property  allows  the  construction  of  grids  over 
contoured  lines  or  surfaces  by  propagating  in  space  the  initial  information  provided  by 
these  lines  or  surfaces.  The  grid  generation  method  must  be  carefully  chosen  to  facilitate 
the  type  of  grid  desired.  Each  type  of  grid  generation  method  will  require  some  addi¬ 
tional  data  to  allow  for  a  suitable  solution.  To  some  extent  this  additional  data  may 
determine  which  type  of  grid  method  should  be  utilized.  If  x  and  y  are  spatial  coordi¬ 
nates  (which  is  true  for  the  present  case)  then  the  additional  data  will  be  the  boundary 
conditions.  If  x  and  y  represent  time  then  the  additional  data  will  represent  initial  con¬ 
ditions  . 

A.  ELLIPTIC  GRID  GENERATION 

Field  grid  generation  obtained  by  the  solution  of  an  elliptic  set  of  equations  requires 
specification  of  each  and  every  boundary  point  of  the  closed  domain  where  the  grid  will 
be  generated.  The  inner  boundary  is  simply  the  body  surface  grid  and  the  outer 
boundary  is  a  user  specified  shape.  The  body  surface  must  be  specified  exactly.  How¬ 
ever,  there  is  some  flexibility  in  choosing  the  shape  of  the  outer  boundary.  Another  re¬ 
quirement  of  this  method  is  that  the  curvilinear  coordinates  must  be  constant  or 
monotonic  on  the  boundaries.  If  any  extrema  of  the  curvilinear  coordinates  exist  in  the 
interior  of  the  physical  region  then  overlapping  of  the  grid  lines  will  occur.  When  using 
an  elliptic  method,  initial  boundary  slope  discontinuities  are  not  propagated  into  the 
field.  This  feature  of  elliptic  grid  generators  tends  to  make  the  grid  very  smooth.  The 
large  computational  time  requirement  for  the  solution  of  the  elliptic  system  of  P.D.E.s 
can  be  a  disadvantage.  The  simplest  form  of  an  elliptic  P.D.E.  is  the  Laplace  equation; 

VY  =  0  (5) 

where  /=  1,2  for  two-dimensional  grid  generation  and  i  =  1,2,3  for  three-dimensional 
grid  generation.  The  effect  of  the  Laplace  operator  is  that  a  very  smooth  grid  is 
produced  which  becomes  equally  spaced  away  from  the  boundary.  The  Laplacian  also 
guarantees  one  to  one  mapping  of  the  coordinate  system.  This  method  will  have  the 
effect  of  making  the  grid  lines  more  closely  spaced  over  concave  boundaries  and  sparse 
over  convex  boundaries.[Ref.  4  :  pp.  18S-22S][Ref.  25:  pp.  503-510] 

Another  approach  to  generate  the  field  grid  is  to  solve  a  Poisson  system  of 
equations.  This  system  has  the  following  general  form; 


35 


vY  =  p‘. 


(6) 


The  forcing  term  P'  can  be  used  to  control  the  spacing  and  orientation  of  the  grid  lines. 
This  control  can  be  extended  to  move  the  intersection  slope  of  the  grid  line  with  the 
boundary.  When  P'  -*  0  the  grid  lines  tend  to  become  equally  spaced,  i.e.  to  approach 
a  grid  obtained  from  the  solution  of  Laplace's  equation.  The  forcing  term  P1  can  also 
be  used  to  enhance  grid  orthogonality.  Orthogonality  of  the  grid  lines  close  to  the  body 
surface  grid  does  not  occur  normally  in  elliptical  solutions.  The  main  advantage  of  a 
Poisson  type  grid  generator  is  that  orthogonality  control  of  the  grid  lines  can  be  main¬ 
tained  at  the  expense  of  complex,  lengthy  and  expensive  calculations.  There  are  several 
elliptic  grid  generators  in  use  today  that  maintain  grid  orthogonality.! Ref.  4:  pp. 
193-236.]  [Ref.  26] 

B.  HYPERBOLIC  GRID 

The  hyperbolic  method  involves  marching  in  space  in  a  time-like  fashion  of  the 
boundary  information,  i.e.,  a  surface  grid.  This  method  is  suitable  for  external  flow 
problems  where  the  exact  location  and  shape  of  the  outer  boundary  is  of  no  vital  im¬ 
portance.  One  major  advantage  is  that  computationally  this  method  is  efficient;  in  ad¬ 
dition.  orthogonality  of  the  field  grid  is  preserved.  Hyperbolic  methods  are  usually  one 
to  two  orders  of  magnitude  faster  than  the  elliptic  methods  because  of  their  noniterative 
time-like  marching  nature.  Control  of  the  grid  lines  is  somewhat  restrictive  but  specifi¬ 
cation  of  the  cell  volume  can  result  in  the  avoidance  of  overlapping  grid  lines,  especially 
in  concave  areas.  Overlapping  of  the  grid  lines  is  not  allowed  because  singularities  are 
propagated  into  the  field,  so  great  care  must  be  taken  to  avoid  these  when  constructing 
the  surface  grid.  Because  the  characteristics  of  the  hyperbolic  method  include 
orthogonality  preservation  and  computational  efficiency,  a  hyperbolic  grid  generation 
method  was  selected  as  opposed  to  an  elliptic  grid  generation  method.  Great  care  was 
taken  in  the  generation  of  the  surface  grid  to  remove  any  singularities  such  as  sharp 
corners,  because  rapid  transitions  on  the  surface  geometry  usually  produce  intersecting 
grid  lines  of  the  field  grid.  The  wing  configuration  did  not  contain  any  severe  concave 
surfaces  which  would  cause  intersecting  of  the  grid  lines.  Two  different  grid  topologies 
were  examined  for  the  grid  generation  of  the  field  grid  over  a  double-delta  wing,  i.e., 
cylindrical  and  spherical.  First  the  cylindrical  grid  generation  procedure  is  presented. 
[Ref.  4:  pp.  272-276][Ref.  25:  p.  503] 
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1.  Cylindrical  Grid  Generation 

The  cylindrical  type  grid  produces  an  H-0  configuration  shown  in  Figure  (23). 
The  grid  in  this  figure  is  called  H-0  because  when  the  wing  is  viewed  from  above  or  from 
the  side  the  erid  appears  to  have  a  "H"  shape.  When  the  grid  is  viewed  from  a  nose-on 
viewpoint  (i.e.  the  yz-plane),  then  the  grid  has  an  "O"  shape,  see  Figure  (24).  It  is  usu¬ 
ally  easier  to  generate  a  cylindrical  grid  than  generating  a  spherical  grid  because  the 
three-dimensional  cylindrical  grid  is  an  assembly  of  planar  two-dimensional  sections. 
This  combination  of  two-dimensional  planes  starts  at  the  nose  and  progresses  aft  to¬ 
wards  the  trailing  edge,  see  Figure  (25)  and  Figure  (26).  Various  methods  may  be  used 
to  generate  these  two-dimensional  grids  on  the  planar  cross-sections.  Here,  a  two- 
dimensional  hyperbolic  grid  generator  was  used  to  generate  the  plane  O-type  grids  at 
various  locations  along  the  axial  direction.  The  cylindrical  grid  has  a  singular  point  at 
the  apex,  Figure  (25),  where  special  care  must  be  taken  during  numerical  solution  for  the 
computation  of  the  transformation  metrics  and  the  application  of  boundary  conditions. 
Along  the  singular  line  starting  from  the  apex  and  extending  upstream  to  the  beginning 
of  the  domain  all  the  grid  points  collapse  on  to  a  single  point. 

Before  the  surface  grid  data  points  generated  for  a  generic  surface  can  be  used 
for  a  cylindrical  grid,  a  modification  was  required  at  the  nose  of  the  grid.  As  can  be  seen 
in  Figure  (7)  a  singularity  is  present  at  the  very  first  point  at  the  apex.  To  avoid  the 
problems  that  this  point  can  cause,  the  singular  point  at  the  apex  was  omitted.  The  field 
grid  was  then  generated  to  have  an  annular  type  appearance  in  the  yz-plane  cross- 
section.  see  Figure  (29).  Because  the  radius  is  so  small,  the  effects  due  to  the  inaccuracy 
of  its  surface  definition  are  negligible. 

The  process  of  improving  the  quality  of  the  field  grid  was  completed  using  the 
computer  graphics  program  PLOT3D  [Ref.  27].  This  graphics  package  was  designed  to 
facilitate  visualization  of  the  field  and  surface  grids  and  flow  fields  of  computed  and  ex¬ 
perimental  results.  This  same  program  was  utilized  to  correct  the  surface  grid  during  its 
development  process.  Because  of  the  early  concern  regarding  the  surface  grid  and 
knowledge  of  hyperbolic  grid  shortcomings,  only  minor  adjustments  were  required  in  the 
surface  grid.  These  adjustments  required  a  redistribution  of  surface  grid  points  in  the 
spanwise  direction. 

The  initial  surface  grid  had  a  linear  distribution  that  was  deemed  inadequate  just 
by  visual  inspection  and  from  physical  considerations  of  the  flow  field  character  at  the 
leading  edge  region.  A  quadratic  distribution  was  subsequently  used  which  appeared  to 
yield  a  better  distribution  of  the  grid  points.  Not  until  the  field  grid  was  generated  was 
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it  clear  that  this  distribution  was  still  insufficient.  The  necessary  changes  in  the  field  grid 
were  to  smooth  out  the  distribution  near  the  leading  edges  and  to  increase  the  grid 
density  where  vortices  were  most  likely  to  occur.  Finally,  a  sinusoidal  distribution  along 
the  spanwise  direction  gave  the  best  results.  Examples  of  various  cross-sections  of  the 
surface  grid  distribution  for  the  wing  can  be  seen  in  Figures  (29),  (31),  (34)  and  (37). 
At  the  trailing  edge  the  same  reasoning  as  for  the  nose  was  carried  out  and  the  grid  was 
not  collapsed  onto  a  single  line  but  instead  retained  a  finite  thickness.  As  a  result  a  very 
small  but  finite  thickness  wake  was  generated. 

The  field  grid  was  completed  by  extending  the  grid  from  the  wing  apex  to  a  lo¬ 
cation  2.0  -  3.0  chord  lengths  upstream,  where  the  freestream  conditions  can  be  applied. 
The  reason  for  this  is  that  the  flow  field  is  affected  upstream  by  the  presence  of  the  wing. 
This  addition  to  the  grid  was  completed  by  repeating  the  very  first  annular  shaped  grid 
and  by  changing  only  the  x-location.  The  yz-plane  remained  the  same  and  it  wras  re¬ 
peated  as  the  x-distribution  gradually  increased  its  Ax  spacing  until  it  reached  a  user 
specified  x-location  upstream.  The  short  program  for  this  additional  grid  can  be  seen  in 
Appendix  F  and  examples  of  the  entire  field  grid  (cylindrical  type)  can  be  seen  in  Figure 
(24)  through  Figure  (41),  Appendix  B.  The  completed  cylindrical  field  grid  dimensions 
are  110x240x6S. 

The  purpose  of  the  grid  generation  is  to  facilitate  a  numerical  solution  to  a 
system  of  P.D.E.s,  and  for  accurate  solutions  to  occur  certain  areas  require  special 
treatment.  A  problem  that  occurs  in  the  computation  of  derivatives  in  the  apex  region 
is  that  differences  are  taken  between  points  that  may  have  different  flow  characteristics; 
i.e.  one  point  may  be  in  the  boundary  layer  and  the  next  point  may  be  outside  of  this 
regime.  This  is  the  reason  why  clustering  of  grid  points  near  the  apex  is  required  to 
avoid  as  much  as  possible  these  inaccuracies.  Clustering  of  the  grid  points  normal  to  the 
surface  of  the  wing  is  naturally  required  to  resolve  the  velocity  gradients  in  the  viscous 
boundary  layer.  For  the  most  part  though,  the  grid  is  aligned  with  the  main  flow  direc¬ 
tion  and  a  numerical  scheme  that  uses  upwinding  can  produce  accurate  results.! Ref.  20] 
2.  Spherical  Grid 

This  particular  method  of  grid  generation  produces  a  C-0  type  grid  as  can  be 
seen  in  Figure  (24).  This  method  is  more  complex  than  the  cylindrical  one  because  the 
entire  three-dimensional  grid  must  be  generated  simultaneously.  The  spherical  grid 
topology  has  one  singularity  that  is  located  at  the  apex  and  propagates  upstream,  see 
Figure  (42)  and  Figure  (43).  The  spherical  grid  is  also  aligned  with  the  main  flow; 
therefore,  an  upwinding  scheme  can  be  used  for  flow  field  solutions.  One  disadvantage 
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of  the  spherical  grid  topology  is  that  the  visualization  of  the  flow  is  a  little  more  difficult 
than  for  the  cylindrical  grid  because  the  {  =  constant  grid  surfaces  are  not  planes  but 
three-dimensional  surfaces.  The  flow  field  solution  tends  to  converge  faster  on  a  spher¬ 
ical  grid  than  in  the  cylindrical  grid. 

The  hyperbolic  grid  generator  for  this  work  utilized  a  cell  volume  hyperbolic  grid 
generation  scheme.  A  coordinate  transformation  to  the  computational  domain  (£ ,  tj,  £) 
was  performed  where  the  body  surface  was  the  boundary  condition  t,{x,y,  z)  =  0  .  The 
field  grid  is  obtained  by  the  solution  of  the  following  nonlinear  system  of  P.D.E.s; 

Vc+V,(  + Vc =  0 
•Vi  +Vc  +  Vi =  0 

-Wi +  -w* +  Vi2;  -  -Wr?  -  Vi2i  -  W; =  A  v- 

where  the  initial  condition  at  £  =  0  are  the  x,y  ,  and  z  coordinates  of  the  body  surface 
[Ref.  20].  The  first  two  equations  are  the  relations  that  preserve  orthogonality  with  re¬ 
spect  to  an  outward  normal  vector  t  .  The  third  equation  is  a  user  specified  volume 
parameter  that  controls  the  cell  size  and  normal  spacing  of  the  grid  points.  The  grid  is 
generated  by  "marching"  in  the  £  direction  and  the  system  of  P.D.E.s  is  solved  by  an 
approximate-factorization,  noniterative,  implicit  finite  difference  scheme.  Even  though 
grid  orthogonality  and  smoothness  are  maintained  the  quality  of  the  field  grid  is  quite 
sensitive  to  the  quality  of  the  surface  grid.  Control  of  grid  clustering  along  the  normal 
to  the  surface  direction  is  provided,  but  there  is  no  accurate  control  in  the  location  of 
the  outer  boundary  due  to  the  marching  type  solution.  The  outer  boundary  for  the 
purposes  of  the  present  work  is  not  crucial  as  long  as  it  extends  2.0  -  2.5  root  chord 
lengths  away  from  the  body. 

The  three-dimensional  hyperbolic  grid  generation  is  very  sensitive  to  the  surface 
grid  definition,  since  the  surface  grid  distribution  is  propagated  in  space  to  generate  the 
three-dimensional  mesh.  This  sensitivity  to  the  initial  conditions  is  the  reason  why  the 
spherical  grid  generation  presented  more  difficulties  than  the  cylindrical  grid  generation. 
The  apex  singularity  in  conjunction  with  the  sharp  angles  created  most  of  the  problems 
in  the  grid  generation  as  far  as  preservation  of  orthogonality  is  concerned.  Because  of 
this  point  singularity,  a  blunting  of  the  nose  was  first  attempted.  This  nose  region  also 
resulted  in  a  reduction  of  grid  lines  in  the  x-direction  to  1 10.  For  the  C-0  configuration 
the  first  grid  line  extends  upstream  and  therefore  it  is  not  necessary  to  add  axial  grid  lines 
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as  in  the  cylindrical  grid  case.  For  effective  Field  grid  generation,  distributions  in  all  three 
directions  had  to  be  adjusted  near  the  nose.  After  many  iterations  a  solution  was  gen¬ 
erated  that  required  only  minor  adjustments.  These  adjustments  were  made  by  writing 
a  small  program  that  would  algebraically  adjust  the  x-distribution  at  the  nose  singularity. 
Algebraically  fixed  were  the  second  and  third  grid  planes  in  the  x-direction.  This  tech¬ 
nique  is  sometimes  the  only  option  available  when  the  grid  requires  minor  adjustments 
and  the  changing  of  usual  parameters  produces  negative  results.  Although  a  suitable 
grid  was  constructed,  an  alternative  method  of  allowing  the  apex  to  collapse  to  a  sharp 
point  was  attempted  in  hopes  of  an  even  more  suitable  field  grid.  Also  at  this  time  sol¬ 
utions  were  being  generated  for  the  cylindrical  grid  and  it  was  observed  that  the  field  grid 
required  a  higher  grid  density  in  the  wing  area  to  facilitate  better  definition  of  the  vortex 
that  develops  over  it.  For  this  reason  the  number  of  grid  lines  in  the  axial  direction  was 
increased  to  160.  By  allowing  the  apex  to  converge  linearly  to  a  point  resulted  in  a 
three-dimensional  mesh  that  had  been  unsurpassed  up  to  this  point.  The  disadvantage 
of  allowing  the  strake  to  linearly  collapse  to  a  point  was  the  original  y-direction 
sinusoidal  distribution  (0°  to  45°)  of  the  surface  grid  was  now  insufficient.  A  solution 
to  this  problem  was  to  allow  the  y-direction  distribution  near  the  apex  to  have  a  45°  to 
60°  distribution.  This  distribution  results  in  a  spacing  of  the  grid  lines  that  i:  nearly 
linear.  Then  the  distribution  was  allowed  to  change  incrementally  as  the  x-location 
changed  to  achieve  a  45°  to  90°  sinusoidal  distribution  at  the  strake  and  wing  junction. 
This  change  was  only  required  in  the  y-direction.  The  conscious  decision  to  allow  a 
singularity  at  the  apex  did  not  affect  the  quality  of  the  field  grid  because  the  topology 
of  the  spherical  C-O  type  grid  requires  that  this  area  collapse  to  a  singular  line. 

Other  areas  of  the  wing  did  not  require  adjustments  because  the  surface  grid  did 
not  propagate  any  problems  into  the  field  grid.  Cross-sections  of  the  grid  can  be  seen 
in  Figure  (45)  through  Figure  (56).  The  only  other  adjustment  made  was  to  stretch  the 
entire  grid  to  a  suitable  distance  from  the  wing.  This  also  was  done  with  a  small  pro¬ 
gram  that  operated  on  the  data  file  that  was  generated  from  the  hyperbolic  grid  genera¬ 
tor.  While  the  final  grid  appears  more  uniform  throughout  the  entire  field,  slight 
deviations  in  the  orthogonality  to  the  surface  did  occur.  The  final  grid  dimensions  for  the 
spherical  field  grid  were  160x240x68  and  can  be  seen  in  Figure  (42)  through  Figure  (56), 
Appendix  C. 
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C.  PARABOLIC  GRID  GENERATION 

The  last  grid  generation  technique  based  on  the  solution  of  P.D.E.s  is  the  parabolic 
grid  generation  method.  The  parabolic  grid  generation  techniques  may  be  constructed 
by  modifying  elliptic  methods  and  hence  carries  various  advantages  of  the  method.  The 
most  popular  modification  is  the  elimination  of  the  second  derivatives.  The  solution  is 
generated  by  marching  out  in  one  direction  like  in  the  hyperbolic  method,  but  the 
marching  is  influenced  somewhat  by  the  other  boundary.  Control  functions  can  be  used 
to  enhance  orthogonality,  which  would  not  occur  normally.  Because  of  the  effect  of  the 
other  boundaries  these  methods  tend  to  have  more  smoothing  effects  than  a  hyperbolic 
method.  The  parabolic  method  has  the  characteristics  that  are  present  in  both  elliptic 
and  hyperbolic  grid  generation  methods.  The  complexity  tends  to  be  less  than  the  el¬ 
liptical  method  and  hence  is  faster.[Ref.  4  :  pp.  277-278] 
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V.  RESULTS  AND  DISCUSSIONS 


Modem  fighter  aircraft  designs  take  advantage  of  the  strakes  to  improve 
controllability  and  enhance  lift  capabilities  at  high  angles  of  attack.  Existing  aircraft  are 
currently  structurally  modified  by  adding  leading  edge  extensions  to  wings  in  order  to 
improve  the  fuselage  flow'  field  characteristics  which  would  eventually  lead  to  improved 
lift  and  maneuverability  characteristics  [Ref.  28]  .  Complex  fluid  dynamic  phenomena 
are  associated  with  the  high  angle  of  attack  flow'  over  the  forebody,  the  strake  and  the 
wings.  The  flow  separates  to  form  vortices  which  provide  nonlinear  lift.  At  high  angles 
of  attack  the  forebody,  the  strake  and  the  wing  vortices  interact  with  each  other  and  as 
a  result  self-excited  unsteady  flow  may  be  triggered.  When  the  angle  of  attack  is  further 
increased  vortex  breakdown  occurs  which  will  enhance  flow-  unsteadiness  and  may  result 
in  a  loss  of  controllability  and  other  undesirable  effects  such  as  wing  rock  or  tail  buffet. 
Many  of  these  interesting  flow  phenomena  can  be  observed  for  the  flow  over  the 
strakc-delta  wing  configuration  model  for  which  the  grid  was  generated. 

In  this  chapter  a  survey  of  the  grid  generation  will  be  done  and  the  results  of  the 
numerical  solution  showing  the  characteristics  of  the  flow  field  will  be  presented.  Dis¬ 
cussions  on  vortex  formation,  interaction  and  breakdown  will  be  made  for  the  various 
angles  of  attack. 

A.  GRID  GENERATION 

The  generation  of  the  surface  and  the  field  grid  is  a  prerequisite  for  the  numerical 
solution.  The  grid  generation  can  be  a  very  time  consuming  process.  However,  the  grid 
quality  will  determine  the  accuracy  of  the  numerical  solution.  The  double-delta  wing 
analyzed  here  is  a  simple  model  of  a  modern  fighter  aircraft  planform.  The  surface  de¬ 
finition  can  be  done  entirely  using  linear  relationships.  This  allowed  to  use  relatively 
simple  algebraic  and  geometric  relationships  to  generate  a  surface  grid.  Even  with  these 
simplifications  the  time  expended  on  creating  a  surface  grid  and  two  types  of  field  grids 
was  quite  long.  The  amount  of  effort  that  is  expended  on  grid  generation  for  an  actual 
aircraft  configuration  would  very  likely  be  more  than  one  year. 

For  symmetric  bodies  it  is  sufficient  to  generate  half  the  surface  grid.  The  gener¬ 
ation  of  the  surface  grid  was  simplified  by  dividing  the  wing  into  similar  sections  and 
writing  computer  programs  specific  for  the  cross-section.  It  was  found  to  be  simpler  to 
progress  from  the  nose  to  the  tail  by  a  A.r  increment  and  compute  the  points  (in  the 
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yz-plane)  that  represent  the  wing  cross-section  at  that  particular  x-location.  The  major 
concern  during  this  phase  was  to  write  a  code  that  would  permit  a  variable  distribution 
of  the  grid  lines.  It  was  also  equally  important  that  ..he  distribution  of  grid  lines  must 
transition  as  evenly  as  possible  between  fine  distributions  and  coarser  distributions.  Fine 
distributions  occur  at  locations  where  large  numbers  of  grid  points  are  required.  Coarse 
distributions  were  chosen  at  locations  where  small  grid  density  suffices  to  capture  the 
physics  of  the  flow.  While  a  finished  surface  grid  may  appear  to  have  smooth  distrib¬ 
utions  and  be  very  uniform,  only  the  subsequent  field  grid  generation  revealed  whether 
the  chosen  distributions  were  adequate.  Therefore,  great  emphasis  should  be  placed  on 
surface  grid  versatility  in  the  early  stages. 

The  distribution  of  the  grid  lines  and  the  number  of  grid  lines  in  the  x-direction  was 
changed  during  the  interactive  process  of  improving  the  surface  grid  to  produce  a  suit¬ 
able  field  grid.  In  the  case  of  the  cylindrical  grid  generation,  this  trial  and  error  process 
was  relatively  short  because  the  field  grid  is  composed  of  two-dimensional  grids.  On  the 
other  hand,  for  spherical  grid  generation  both  axial  and  circumferential  surface  grid  dis¬ 
tributions  were  much  more  sensitive  because  they  must  be  suitable  for  the  generation 
of  a  three-dimensional  mesh.  During  the  process  of  generating  both  the  cylindrical  and 
spherical  type  grids,  various  small  programs  were  written  to  refine  and  improve  a  par¬ 
ticular  grid  that  was  generated.  Constructing  a  grid  requires  a  knowledge  and  familiarity 
with  grid  generation  codes  and  some  prior  knowledge  of  the  flow  characteristics.  For 
example,  the  cylindrical  grid  required  the  deletion  of  the  first  points  of  the  original  sur¬ 
face  grid  in  order  to  eliminate  the  singularity  at  the  apex.  Another  program  was  then 
written  to  extend  the  grid  upstream  to  where  conditions  of  the  freestream  where  ex¬ 
pected.  These  programs  can  be  found  in  Appendix  F.  The  generation  of  the  spherical 
grid  required  the  writing  of  additional  small  programs  in  an  attempt  to  algebraically 
modify  regions  of  the  field  grid  that  had  small  discontinuities,  see  Appendix  F.  Diligence 
in  changing  the  distributions  and  alteration  of  the  nose  region  resulted  in  a  smooth 
spherical  field  grid.  It  is  believed  that  despite  the  larger  amount  of  time  spent  for  the 
generation  of  the  splverica!  grid,  a  better  quality  grid  compared  with  the  cylindrical  grid 
was  obtained.  However,  because  of  a  time  constraint  this  could  not  be  verified  by  ob¬ 
taining  a  solution  on  the  spherical  grid.  The  cylindrical  grid  was  used  for  the  numerical 
implementation  because  a  flow  solution  was  desired  for  the  presentation  of  this  work. 

A  visual  comparison  of  the  two  completed  field  grids  reinforces  the  opinion  that  a 
spherical  grid  topology  is  more  suitable  for  the  subsequent  numerical  solution.  It  is 
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emphasized  however,  that  spherical  grid  generation  is  very  sensitive  to  the  surface  grid 
distributions  and  smoothness  and  requires  larger  computing  times. 

B.  FLOW  FIELD  CHARACTERISTICS 

The  main  characteristic  of  the  flow  over  delta  wings  is  the  presence  of  the  leading 
edge  vortex.  The  nonlinear  induced  lift  by  the  leading  edge  vortices  has  been  actively 
investigated  in  recent  years.  Vortical  flow  is  an  advantageous  lift  generation  mechanism 
that  can  be  utilized  successfully  at  medium  to  high  speeds.  A  description  of  the 
leeward-side  flow  field  characteristics  will  be  presented  and  the  flow  field  structure  over 
a  double-delta  wing  at  various  angles  of  attack  will  be  analyzed.  Available  experimental 
results  will  be  used  to  validate  the  computed  results. 

1.  Vortex  Characteristics 

The  leading  edge  vortices  result  from  the  roll-up  of  the  shear  layer  that  is  shed 
from  the  leading  edge.  At  moderate  to  high  angles  of  attack  the  wingward  and  leeward 
side  flow  merges  to  form  a  shear  layer  that  rolls  up  and  forms  a  rotational  vortex  core. 
In  the  case  of  the  double-delta  wing  a  system  of  two  primary  vortices  is  formed.  The 
first  is  along  the  strake  and  the  second  along  the  leading  edge  of  the  wing  section.  The 
primary  strake  vortex  reattaches  itself  to  the  centerline  of  the  planform.  The  vortex 
strength  is  increased  by  a  continuous  feeding  of  vorticitv  from  the  shear  layers  of  the 
leading  edge.  Surface  pressure  suction  peaks  are  produced  at  a  location  below  the  po¬ 
sition  of  the  vortex  core.  Because  the  vortex  core  exhibits  large  gradients  of  vorticitv 
and  circumferential  velocity,  large  viscosity  effects  are  expected.  As  the  vortex  strength 
increases  downstream  so  do  the  lateral  velocities  that  are  near  the  surface.  Coincidental 
with  large  velocities  is  a  decrease  in  pressure.  A  secondary  separation  and  formation  of 
the  secondary  vortex  is  the  result  of  these  large  lateral  velocities  and  the  associated  ad¬ 
verse  pressure  gradients.  If  the  secondary  vortex  is  strong  enough,  a  tertiary  vortex  can 
form  under  the  secondary  vortex  by  the  same  mechanisms  [Ref.  17).  A  schematic 
showing  the  leeward  side  vortex  system  and  sense  of  rotation  of  these  vortices  is  shown 
in  Figure  (6).  Separated  flow  of  the  secondary  vortex  system  reattaches  again  on  the 
wing  leeward  surface.  This  separation  and  reattachment  process  can  be  best  visualized 
with  the  use  of  surface  oil  flow  patterns,  see  Figure  (57).  This  visualization  method  has 
effectively  shown  the  separation  and  reattachment  of  both  primary,  secondary  and  ter¬ 
tiary  vortices. 

The  strength  of  the  leading  edge  vortex  increases  downstream  and  as  the  angle 
of  attack  is  increased.  The  pressure  gradient  in  the  direction  of  a  vortex  core  accelerate 
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Figure  6.  Primary,  Secondary  and  Tertiary  Vortices 

the  fluid  particles  until  a  critical  angle  of  attack  is  reached.  At  this  angle  the  organized 
vortex  core  suddenly  breaks  down  due  to  the  adverse  pressure  gradient  at  the  trailing 
edge.  .  This  sudden  transition  is  more  commonly  known  as  vortex  burst  or  vortex 
breakdown.  Vortex  burst  is  a  flow  phenomenon  that  needs  to  be  understood  because  a 
loss  of  the  suction  peak  will  occur  and  this  change  of  the  induced  lift  can  result  in  un¬ 
desirable  effects  on  the  aircraft.  Work  in  this  area  by  Sarpkaya,  Thomas,  Kjelgaard  and 
Sellers.  Ekaterinaris,  Hawk.  Barnett  and  O'Niel,  Kegelman  and  Roos  and  many  others 
has  been  conducted  in  recent  years  (Ref.  12,29,30,20,31,32]. 
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Vortex  burst  is  a  transition  from  a  jet-like  spiraling  flow  to  a  wake-like  flow. 
Adverse  pressure  gradient  and  swirl  angle  of  the  flow  contribute  to  this  phenomenon. 
Swirl  angle  in  defined  as; 


4>  =  tan-’(-^) 


where  u  is  the  axial  velocity  component  and  ve  is  the  azimuthial  velocity  component. 
Vortex  breakdown  will  usually  occur  when  the  swirl  angle  exceeds  a  critical  value  of 
approximately  40°  .  The  swirl  angle  and  the  adverse  pressure  gradient  determine  the 
type  of  burst  for  a  cylindrical  vortex,  i.e.,  bubble  breakdown,  spiral  breakdown  or  double 
helix  breakdown  [Ref.  12].  Normally,  as  the  angle  of  attack  is  increased  the  burst  lo¬ 
cation  will  move  upstream.  If  the  angle  of  attack  is  increased  even  more,  wake  type  of 
flow  behind  a  bluff  body  will  be  encountered.  Most  of  the  initial  investigations  of  vortex 
generation,  induced  lift  and  vortex  breakdown  was  completed  using  a  single-delta  wing 
configuration.  In  this  investigation  the  complex  flow  field  that  results  due  to  the  pres¬ 
ence  of  multiple  vortices  is  even  further  complicated  for  the  case  of  the  double-delta 
wing.  This  is  due  to  the  interaction  of  the  strake  vortex  with  the  wing  vortex  and  with 
the  surface  of  the  wing.  The  characteristics  are  shown  for  various  angles  of  attack  and 
results  are  discussed  below  in  more  detail. 

2.  Double-Delta  Wing  Flow  Characteristics 
a.  Angle  of  Attack  -  10° 

The  computed  surface  flow  pattern  at  a  =  10°  shown  in  Figure  (57)  does  not 
indicate  tertiary  separation  on  the  strake,  while  secondary  and  tertiary  separation  are 
shown  on  the  wing.  The  leeward  side  flow  characteristics  of  the  double-delta  wing  show 
at  10.0°  angle  of  attack  are  shown  in  Figure  (58)  and  Figure  (59).  Two  primary  vortices 
are  formed,  the  first  is  formed  by  the  sharp  leading  edge  of  the  strake  and  the  second 
by  the  leading  edge  of  the  wing.  Both  of  these  vortices  are  continually  fed  by  vorticity 
as  they  progress  downstream  which  increases  their  strength.  The  sense  of  rotation  for 
the  primary  strake  vortex  and  primary  wing  vortex  is  the  same  as  can  be  seen  in  the  ve¬ 
locity  vector  diagrams  in  Figure  (60)  through  Figure  (62).  Also  clearly  shown  in  these 
figures,  are  the  primary  and  secondary  vortices  having  opposite  swirling  directions. 
Depending  on  the  angle  of  attack,  both  primary  vortices  may  swirl  around  each  other, 
but  here  the  vortices  remain  separated.  On  the  other  hand,  the  wing  tip  vortex  and  the 
wing  vortex  do  eventually  merge  as  can  be  seen  in  Figure  (58). 
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At  10°  angle  of  attack  the  primary  strake  vortex  separates  and  reattaches 
at  the  centerline  of  the  wing,  see  Figure  (57).  It  can  also  be  seen  that  a  secondary  sep¬ 
aration  occurs  near  the  leading  edge  of  the  strake.  This  secondary'  vortex  also  reattaches 
itself,  but  it  is  not  strong  enough  at  10°  to  generate  a  tertiary  vortex.  The  wing  vortex 
and  reattachment  of  the  primary  and  secondary  vortex  can  also  be  seen  in  Figure  (57). 

b.  Angle  of  Attack  -  19° 

At  this  angle  of  attack  the  primary  strake  vortex  is  much  stronger.  This  can 
be  deduced  by  the  presence  of  a  tertiary  vortex  as  seen  in  Figure  (63).  This  increase  of 
vortex  strength  would  produce  an  increase  of  the  induced  lift.  The  wing  vortex  that  de¬ 
velops  is  continually  fed  by  two  sources.  One  source  is  the  shear  layer  that  is  connected 
to  the  wing  leading  edge  and  the  other  is  the  shear  layer  associated  with  the  primary 
strake  vortex.  This  relinquishment  of  vorticity  by  the  primary'  strake  vortex  causes  its 
strength  downstream  of  the  kink  to  remain  constant  or  even  to  reduce.  These  two 
vortices  eventually  merge  close  to  the  trailing  edge  of  the  wing,  see  Figure  (64).  The 
vortex  burst  defines  the  limit  of  vortex  strength  that  can  be  maintained  by  the  flow  field. 
The  burst  appears  to  occur  shortly  after  the  primary  strake  and  primary'  wing  vortices 
merge,  Figure  (64)  and  Figure  (6.5).  It  is  at  this  angle  of  attack,  i.e.  just  before  vortex 
burst  appears,  where  the  maximum  induced  lift  occurs.  As  the  angle  of  attack  is  in¬ 
creased  further  the  strake  vortex  continues  to  get  stronger  but  the  burst  point  moves 
further  upstream.  Close  examination  of  particle  traces,  see  Figure  (64)  and  Figure  (65), 
actually  shows  the  development  of  a  wing-tip  vortex.  This  wing-tip  will  eventually  merge 
with  the  primary  wing  vortex.  The  flow  direction  for  both  the  strake  and  wing  primary 
vortices  is  the  same,  see  Figure  (66)  through  Figure  (68).  These  cross-sectional  views 
of  the  velocity  vectors  also  show  the  counter  rotation  between  primary,  secondary  and 
tertiary  vortices.  At  the  cross-section  over  the  wing,  Figure  (66),  the  two  vortex  cores 
are  distinct.  While  at  the  trailing  edge,  Figure  (68),  the  two  vortices  are  merged. 

c.  Angle  of  Attack  -  22.4° 

As  the  angle  of  attack  continues  to  increase  the  most  prominent  feature  is 
probably  the  change  in  location  of  the  vortex  breakdown.  The  location  of  the  burst  will 
move  upstream  as  the  angle  of  attack  is  increased.  Not  only  does  the  burst  location 
move  upstream,  but  the  strength  of  the  strake  vortex  increases.  Figure  (69)  clearly 
shows  the  development  of  a  tertiary  vortex  which  is  characteristic  for  a  strong  vortex. 
The  breakdown  of  the  vortex  is  readily  apparent  in  the  22.4°  flow  solutions.  In  Figure 
(70)  and  Figure  (71)  the  bursting  of  the  strake  vortex  over  the  wing  can  be  seen.  The 
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rotations  of  the  vortices  that  develop  here  are  the  same  as  at  smaller  angles  of  attack 
and  can  be  verified  by  Figure  (72)  through  Figure  (74). 

3.  Comparison  with  Experimental  Data 

Although  the  purpose  of  Cunningham  and  Boer's  experiment  was  the  investi¬ 
gation  of  unsteady  phenomena,  the  report  also  presents  steady  state  data.  The  devel¬ 
opment  and  location  of  the  vortices  are  in  qualitative  agreement  with  Cunningham  and 
Boer's  flow  visualization  results  (Ref.  18].  These  authors  also  present  steady  pressure 
measurements  for  five  different  angles  of  attack.  Unfortunately,  insufficient  time  was 
available  toward  the  completion  of  this  investigation  to  attempt  a  detailed  comparison 
of  the  present  computational  results  with  the  pressure  data.  The  only  comparison  made 
was  for  the  pressures  calculated  at  a.  =  19°.  Similarly,  the  double-delta  wing  studied  by- 
Krause  and  Liu  [Ref.  17]  and  the  experimental  data  of  Brennenstuhl  cited  therein  has 
not  yet  been  used  for  comparison  purposes. 

A  comparison  of  tne  pressure  coefficient  and  spanwise  location  to  experimental 
data  can  be  seen  in  Figure  (75)  through  Figure  (77).  Three  axial  locations  on  the  wing 
were  selected,  specifically,  xc  =  0.40.  0.66,  0.98.  The  locations  were  selected  to  inves¬ 
tigate  representative  sections  of  the  strake,  the  wing  and  the  trailing  edge.  It  can  clearly 
be  seen  that  at  x  c  =  0.40  there  develop  two  suction  peaks  due  to  the  primary  and  sec¬ 
ondary  vortex.  The  location  of  these  peaks  differs  from  the  experiment,  but  the  trend 
is  well  represented.  Since  the  calculations  are  believed  to  be  not  yet  fully  converged,  it 
is  expected  that  more  computational  time  will  yield  even  more  accurate  results.  It  must 
be  pointed  out  that  these  calculations  are  for  vortical  separated  flow  and  the  accurate 
capture  of  these  characteristics  is  very  difficult.  In  Figure  (76),  a  cross-section  of  the 
wing  shows  the  two  suction  peaks  due  to  the  primary  strake  and  primary  wing  vortices 
quite  well.  As  before,  the  calculations  reproduce  the  trends  quite  well,  but  it  is  expected 
that  more  fully  converged  results  will  produce  still  better  agreement.  The  last  compar¬ 
ison  is  made  at  the  trailing  edge  and  is  shown  in  Figure  (77).  Again,  the  trend  is  re¬ 
produced  well:  but,  because  at  this  axial  location  vortex  breakdown  occurs,  it  is  very- 
difficult  to  achieve  exact  results.  This  graph  shows  a  smaller  pressure  coefficient  com¬ 
pared  to  the  other  axial  locations.  This  decrease  is  a  direct  result  of  the  vortex  break¬ 
down. 

All  the  comparisons  of  the  coefficient  of  pressure  depicted  quite  well  the  trends 
associated  with  separated  vortical  flow.  It  is  expected  that  the  locations  of  the  suction 
peaks  will  be  even  more  accurate  after  more  convergence.  The  comparison  made  here 
are  based  on  results  that  required  35-40  hours  of  CPU  time  on  the  Cray-Y.MP  computer. 
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It  is  expected  that  another  ten  hours  of  CPU  time  would  produce  even  more  accurate 
results.  The  expected  total  CPU  time  for  each  angle  of  attack  is  approximately  50  hours. 
Due  to  insufficient  time,  the  fully  converged  results  could  not  be  presented  in  this  paper. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


An  investigation  of  the  flow  characteristics  that  are  created  in  the  fluid  domain  sur¬ 
rounding  a  double-delta  wing  at  various  angles  of  attack  by  a  numerical  approach  was 
presented.  The  numerical  solution  of  the  .Navier-Stokes  equations  requires  the 
discretization  of  the  flow  domain  by  a  smooth  computational  grid.  Accurate  represen¬ 
tation  of  the  flow  physics  by  the  grid  points  will  directly  affect  the  quality  of  the  flow 
field  solution.  Therefore,  the  surface  and  field  grid  density  and  topology  must  be  care¬ 
fully  chosen.  With  the  flow  field  domain  defined  the  numerical  solution  (in  this  case  fi¬ 
nite  difference)  can  be  implemented  to  investigate  the  flow  characteristics  or  compare 
with  experimental  results. 

An  algebraic  method  of  grid  generation  was  selected  for  the  defining  of  the  surface 
grid  and  the  source  code  is  presented  in  Appendix  E.  The  important  precaution  is  that 
when  developing  the  source  code  attempt  to  allow  the  flexibility  of  as  many  parameters 
as  possible.  The  grid  line  distribution  and  number  of  grid  points  in  a  specified  direction 
were  found  to  be  most  important.  Subsequent  generation  of  the  field  grid  will  require 
the  moving  of  grid  lines  to  a  distribution  that  will  produce  a  smooth  and  continuous 
grid. 

Different  types  of  numerical  grid  generation  techniques  are  available,  such  as 
hyperbolic,  elliptic  and  parabolic.  The  advantages,  disadvantages  and  characteristics  of 
these  methods  were  discussed.  The  hyperbolic  grid  generation  technique  was  chosen  and 
two  field  grid  topologies  were  generated,  a  cylindrical  grid  and  a  spherical  grid.  The  cy¬ 
lindrical  grid  was  easier  to  generate,  but  the  spherical  grid  yielded  a  smoother  grid  dis¬ 
tribution  in  space.  This  was  achieved  at  the  expense  of  time  and  computational  effort. 
However,  use  of  the  cylindrical  grid  would  allow  the  generation  of  an  acceptable  flow 
field  solution. 

Once  the  grid  was  completed,  a  finite  difference  algorithm  in  conjunction  with  an 
algebraic  turbulence  model  was  utilized  to  obtain  flow  field  solutions  at 
a  =  10.0°,  19.0°,  22.4°  angles  of  attack.  Investigations  of  vortex  generation,  vortex 
interaction  and  vortex  breakdown  were  conducted.  At  moderate  angles  of  attack  the 
double-delta  wing  configuration  showed  primary  vortices  generated  from  both  the  strake 
and  wing.  These  vortices  produce  nonjinear  vortical  lift  which  can  be  very  beneficial  to 
fighter  type  aircraft  that  operate  at  high  angles  of  attack.  At  approximately  19°  angle 
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of  attack  vortex  bursting  occurred  just  after  the  primary  strake  vortex  and  primary  wing 
vortex  merged  together.  The  formation  of  the  primary  and  secondary  vortices  over  the 
double-delta  wing  compared  favorably  with  the  flow  visualization  data  of  Cunningham 
and  Boer  [Ref.  lSj.  It  is  strongly  recommended,  as  the  next  phase  of  this  investigation, 
to  compare  the  present  numerical  results  with  the  steady  pressure  data  obtained  by  these 
two  authors.  Using  the  source  code  presented  here  as  a  building  block,  future  studies 
could  repeat  the  calculations  for  the  double-delta  wing  studied  by  Krause  and  Liu  [Ref. 
17]  and  then  compare  with  the  experimental  results  of  Brennenstuhl  cited  in  Reference 
17. 

Future  studies  of  this  phenomenon  might  be  to  continue  this  same  analysis  utilizing 
the  spherical  grid  to  determine  if  the  results  are  more  accurate  or  if  computational  time 
is  less,  i.e.  the  solution  field  converges  faster.  Furthermore  with  the  recent  increased 
interest  in  dynamic  stall  phenomenon,  an  analysis  could  be  done  to  compare  the  com¬ 
putational  results  of  a  pitching  straked-wing  to  the  experimental  studies  done  by 
Cunningham  and  Boer.  Because  a  major  portion  of  time  is  spent  in  the  generation  of  a 
field  grid  for  an  analysis  of  this  type,  the  work  load  would  be  reduced  by  using  the  grid 
presented  in  this  thesis. 

The  work  load  involved  in  generating  a  field  grid  is  significantly  increased  when  the 
body  under  investigation  is  an  entire  aircraft.  For  this  reason  the  need  for  a  method  of 
quickly  producing  a  surface  grid  would  expedite  a  numerically  generated  solution  of  a 
flow  field  and  allow  more  time  for  the  improvement  of  numerical  methods. 
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APPENDIX  A.  SURFACE  GRID  FIGURES 


Figure  8.  Double-Delta  Wing  Surface  Grid  (120x240x1) 


Figure  9.  Double-Delta  Wing  -  top  view  (100x240x1) 


Figure  10.  Double-Delta  Wing  and  Wake  -  top  view  (I20.\240\l) 
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Figure  II.  Detail  of  Apex 
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Figure  12.  Grid  Distribution  of  the  Strake  (48x240x1) 
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Figure  14.  Leading  Edge  Rounding  of  the  Strake 
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Figure  15.  Grid  Distribution  of  the  Wing  (26x240\l) 
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figure  16.  Typical  Section  of  the  Wing 
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Figure  18.  Grid  Distribution  of  the  Rectangular  Section  ( I6\240xl) 
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Figure  20.  Leading  Edge  Rounding  of  the  Rectangular  Section 
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Figure  21.  Grid  Distribution  of  the  Wake  (29x121x1) 
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Figure  23.  Cylindrical  (H-O)  Grid  Topology  (130x240x68) 
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Figure  24.  Spherical  (C-O)  Grid  Topology  (160x240x68) 
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APPENDIX  B. 


FIELD  GRID  FIGURES  -  CYLINDRICAL 
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figure  28.  Typical  Cross-section  Upstream  of  the  Apex  -  front  view 


Figure  29.  First  Cross-section  of  the  Strake  -  front  view 
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Figure  31.  Near  Field  Grid  of  the  Strake  -  front  view 


Figure  32.  Leading  Edge  Detail  of  the  Strake  -  front  view 


Figure  33.  Typical  Cross-section  of  the  Wing  -  front  view 
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Figure  37.  Near  Field  Grid  of  the  Rectangular  Section  -  front  view 


Figure  39.  Typical  Cross-section  of  the  Wake  -  front  view 
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Figure  40.  Near  Field  Grid  of  the  Wake  -  front  view 
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Figure  41.  Edge  Detail  of  the  Wake  -  front  view 


APPENDIX  C. 


FIELD  GRID  FIGURES  -  SPHERICAL 
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Figure  42.  Spherical  Grid  Configuration  (160x240x68) 


Figure  43.  Spherical  Grid  Configuration  Detail  (160x240x68) 


Figure  44.  Spherical  Grid  -  side  view  (160x240x68) 
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Figure  46.  Near  Field  Grid  of  the  strake  -  front  view 


Figure  48.  Typical  Cross-section  of  the  Wing  -  front  view 
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Figure  49.  Near  Field  Grid  of  the  Wing  -  front  view 
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Figure  52.  Near  Field  Grid  of  the  Rectangular  Section-  front  view 


Figure  54.  Typical  Cross-section  of  the  Wake  -  front  view 
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Figure  55.  Near  Field  Grid  of  the  Wake  -  front  view 


104 


Figure  56.  Edge  Detail  of  the  Wake  -  front  view 
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Figure  57.  Surface  Flow  Pattern  at  10°  -  M  =  0.22,  Re=  3.8xlOE6,  (70x63x68) 
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Figure  58.  Particle  Traces  at  10°  -  M  =  0.22,  Re  =  3.8\I0E6,  (70x63x68) 
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Figure  59.  Vortex  Location  at  10*  -  M  =  0.22,  Re  =  3J?\I0E6,  (70x63x68) 
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Figure  61.  Wing  Velocity  Vectors  at  10°  -  M  =  0.22,  Re=  3.8\I0E6 


Figure  62.  T.  E.  Velocity  Vectors  at  10°  -  M  =  0.22,  Re  =  3.8xlOE6 
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Figure  63.  Surface  Flow  Pattern  at  I9%-  M  =  0.22,  Re  =  3.8xlOE6,  (70x63x68) 
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Figure  64.  Particle  Traces  at  19°  -  M  =  0.22,  Re  =  3.8.x I0E6,  (70x63x68) 
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Figure  65.  Vortex  Location  at  19°.  -  M  =  0.22,  Re  =  3.8x1 0E6,  (70x63x68) 
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Figure  66.  Strake  Velocity  Vectors  at  19°  -  M  =  0.22,  Re=  3.8.\I0E6 
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Vectors  at  19°  -  M  =  0.22,  Re  =  3.8x1 0E6 
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Figure  68.  T.  E.  Velocity  Vectors  at  19°  -  M  =  0.22,  Re  =  3.8\1 0E6 
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Figure  70.  Particle  Traces  at  22.4"  -  M  =  0.22,  Re=  3.8xlOE6,  (70x63x68) 
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Figure  71.  Vortex  Location  at  22.4°  -  M  =  0.22,  Re  =  3.8x1 0E6,  (70x63x68) 
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Figure  72.  Strake  Velocity  Vectors  at  22.4°  -  M  =  0.22,  Re=  3.8x1 0E6 
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Figure  73.  Wing  Velocity  Vectors  at  22.4°  -  M  =  0.22,  Re  -  3.8xlOE6 


Figure  74.  T.  E.  Velocity  Vectors  at  22.4°  -  M  =  0.22,  Re=  3.8.\I0E6 
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Figure  75.  Surface  Pressure  Coefficient  at  x/c  =  0.40 
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Figure  76.  Surface  Pressure  Coefficient  at  x/c  =  0.66 
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Figure  77.  Surface  Pressure  Coefficient  at  x/c  =  0.9<S 
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nonononnnnonnonnononnn  ooooooooooooooooooooooooooooooo 


THIS  IS  A  PROGRAM  TO  GENERATE  A  SURFACE  GRID  FOR  A  DOUBLE  DELTA 
WING  AIRFOIL  THAT  WILL  BE  ANALYZED  FOR  MY  MASTERS  THESIS. 

CERTAIN  DATA  WILL  BE  REQUIRED  VIA  A  DATA  FILE  AND  INCLUDES  THE 
FOLLOWING  VARIABLES: 

NOSERAD  -  THE  RADIUS  DECSCRIBING  THE  NOSE 

NOSGDX1  -  THE  NUMBER  OF  GRID  POINTS,  IN  THE  X  DIRECTION  DESIRED 
IN  THE  ROUNDED  NOSE 

NOSGDX2  -  THE  NUMBER  OF  GRID  POINTS,  IN  THE  X  DIRECTION  DESIRED 
IN  THE  FIRST  DELTA  WING.  DO  NOT  INCLUDE  IAST  GRID  IN 
THE  NOSE  ROUNDING 

NOSGDY  -  THE  NUMBER  OF  GRID  POINTS,  IN  THE  Y  DIRECTION  DESIRED 
IN  THE  NOSE  BODY 

CRVGDY  -  THE  NUMBER  OF  GRID  POINTS,  IN  THE  Y  DIRECTION  DESIRED 
IN  THE  ROUNDED  LEADING  EDGE 

LEN1  -  THE  LENGTH  FROM  THE  LEADING  EDGE  TO  THE  SECOND  WING 
LEN2  -  THE  HALF  WIDTH  OF  THE  FIRST  DELTA  WING 

LEN3  -  THE  HALF  THICKNESS  OF  THE  FIRST  DELTA  WING 

THE  FOLLOWING  IS  A  DEFINITION  OF  SOME  OF  THE  VARIABLES  USED  IN 
THE  FIRST  PART  OF  THE  PROGRAM: 

DELANG  -  THE  DELTA  ANGLE  USED  TO  GENERATE  GRIDS  IN  THE  Y  DIRECTION 

DISTA  -  THE  DISTANCE  FROM  THE  TIP  OF  THE  AIRFOIL  TO  THE  LAST 

GRID  GENERATED  IN  THE  NOSE  ROUNDING 
DELDIST  -  THE  DELTA  DISTANCE  IN  THE  X  DIRECTION  USED  IN  ROUNDING 
THE  NOSE 

ANG  -  THE  ANGLE  USED  TO  CALCULATE  THE  SPECIFIC  GRIDS  IN  THE 
Y  DIRECTION 

RAD  -  THE  RADIAL  DISTANCE  USED  TO  CALCULATE  GRIDS  IN  THE  Y 
DIRECTION 

XI  -  THE  DISTANCE  FROM  THE  CENTER  OF  THE  NOSE  TO  THE  TRAILING 

EDGE  OF  THE  FIRST  DELTA  WING 

ANG1  -  THE  ANGLE  BETWEEN  THE  CENTERLINE  AND  TRAILING  EDGE 

ANG2  -  THE  ANGLE  BETWEEN  THE  TRAILING  EDGE  AND  LEADING  EDGE 

OF  THE  FIRST  DELTA  WING 

THETA  1  -  THE  ANGLE  FORMED  BY  THE  NOSE  ROUNDING,  IT  WILL  BE 

PERPENDICULAR  TO  THE  LEADING  EDGE  OF  THE  FIRST  DELTA  WING 
DISTB  -  THE  DISTANCE  FROM  THE  TRAILING  EDGE  TO  THE  NOSE  RADIUS 

INUM  -  THE  TOTAL  NUMBER  OF  LINEAR  GRIDS  ON  THE  UPPER  AND  LOWER 

SURFACE 

DISTF  -  THE  DISTANCE  IN  THE  Z  DIRECTION  TRAVERSED  BY  THE  THICKNESS 
OF  THE  LEADING  EDGE 

THETA2  -  THE  ANGLE  FORMED  BY  THE  THICKNESS  OF  THE  FIRST  DELTA  WING 

DISTC  -  THE  HALF  WIDTH  IN  THE  Y  DIRECTION  USED  TO  GENERATE  Y  AND 

Z  DIRECTION  GRIDS 

DISTD  -  THE  HALF  THICKNESS  USED  IN  SAME  CALCULATIONS  AS  DISTC 

THETA 3  -  SAME  TYPE  OF  ANGLE  AS  THETA1  BUT  USED  FOR  Y  AND  Z  GRIDS 

DISTE  -  THE  DISTANCE  IN  THE  Y  DIRECTION  FOR  Y  AND  Z  GRID 

DELY  -  THE  DELTA  Y  DISTANCE  FOR  THE  DISTE  VARIABLE 


DIMENSION  X(170, 150,1) , Y (170, 150 , 1) , Z ( 170 , 150 , 1) 

INTEGER  NOSGDX1 , NOSGDX2 , NOSGDY , CRVGDY , BODGDX1 , BODGDX2 , NOSGDX , 
*  SIDGRD,WAKGPO, CRVGDY 1 

REAL  LEN1 , LEN2 ,  LEN3 ,  LEN4 , LEN5 , LEN6 , NOSERAD, NOSERAD1 
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C 


OPEN  (Li HI T=  10 ,  FILE-'ddwsg.  in' ,  STATUS* 1  OLD' ) 
OPEN ( UNIT” 1 2 , FILE= ' ddwsg . dat ' , STATUS- ' NEW ' ) 


READ ( 10 , *) 
READ (10, * ) 
READ (10,*) 
READ ( 10 , * ) 
READ ( 10 , *) 
READ ( 10 , *) 
READ (10,*) 
READ ( 10 , *) 
READ (10,*) 


NOSERAD , NOSGDX1 , NOSGDX4 , N0SGDX5 , NOSGDY 

CRVGDY , SIDGRD, WAKGRD 

BODGDX1 , BODGDX2 , THETA4 

LEN1 , LEN2 , LEN3 

LEN4 , LEN5 , LEN6 

DISTX1A, DISTX1B, DISTX2 , DISTX3  ,  DISTX4 
CIA , E1A, C1B, E1B, C2A, E2A, C2B, E2B 
C3A,E3A,C3B,E3B 
WAKLEN 


PRINT  * , *  YES1 ' 


THIS  SECTION  DOES  SOME  PRELIMINARY  CALCULATIONS  ON  THE  WING 


PI =2 . 0*ASIN (1.0) 

DELANG=PI/ (NOSGDY-1 ) 

Xl=( (LEN1-N0SERAD)**2. 0+LEH2**2.0)  **0.5 

DUMMY=LEN2/X1 

ANG1=ASIN (DUMMY) 

DUMMY =N0SERAD/X1 
ANG  2  =ACOS ( DUMMY ) 

THETA1=PI-ANG1-ANG2 


THIS  SECTION  DOES  THE  ROUNDING  OF  THE  NOSE  TO  AVOID  SINGULARITIES 


DO  10  IX=1 ,N0SGDX1 

DO  10  IY=1, NOSGDY 

IF (IX. EQ. 1) THEN 

X (IX, IY, 1) -0 . 0 
Y (IX, IY, l)-0. 0 
Z  (IX,  IY ,  1)  «=0 . 0 

ELSE 

X(IX, IY, 1) -NOSERAD-COS ( (THETA1/ (NOSGDX1-1) ) * (IX-1) ) 
*  *NOSERAD 

ANG-PI/2 . 0+ (IY-1) *DELANG 

RAD-NOSERAD*SIN( (THETA1/ (N0SGDX1-1) ) * (IX-1) ) 
Y(IX,IY,1)--1. 0*COS (ANG) *RAD 
Z(IX,IY,1)-SIN (ANG) *RAD 

ENDIF 


10  CONTINUE 

PRINT  *, ' YES2  ' 


THIS  SECTION  DOES  SOME  PRELIMINARY  CALCULATIONS  ON  THE  WING 
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DISTB=LEN1-N0SERAD+C0S (THETA1) *NOSERAD 

INUM-NOSGDY-CRVGDY 

DISTF-LEN  3 -SIN (THETA1 } ‘NOSERAD 

THETA2-ATAN (DISTF/DISTB, 

CUMDIST-0.0 


C 

C  THIS  SECTION  COMPUTES  GRID  POINTS  FOR  THE  NOSE  UP  TO  THE  SECOND  WING 


CRVGDY1-NOSGDY-104 

NOSGDX2=NOSGDX4/2+NOSGDX5/2 

DO  20  IX=l,NOSGDX2 

IF (IX. LT . 4  )  CRVGDY1-CRVGDY1-2 

IF(IX.GE.4 .AND. IX. LT. 6)  CRVGDY1-CRVGDY1-2 

IF ( IX . GE . 6 . AND . IX . LT . N0SGDX2 )  CRVGDY 1-CRVGDY1-2 

I F ( CRVGDY 1 . LT . 7 )  CRVGDY 1-7 

INUMN=NOSGDY-CRVG DY1 

CUMDELY-0. 0 

IF (IX. LE. NOSGDX4/2) THEN 

CALL  STRCH4 ( DISTB , NOSGDX4 , IX, DELDIST, DISTX1A) 
ELSE 


CALL  STRCH4 ( DISTB, NOSGDXS , IX, DELDIST, DISTX1B) 
ENDIF 


ICOUNT=0 

CUMDI ST-CUMDI ST+DELDIST 

DISTC-SIN (THETA1 ) *NOSERAD+ (CUMDIST) /TAN (THETA1 ) 
DISTD-S IN ( THETA 1) *NOSERAD+ ( CUMDI ST) *TAN ( THETA2 ) 
XI- ( (DISTC-NOSERAD) **2 . 0+DISTD**2 . 0) **0 . 5 
DUMMY-DISTD/X1 
ANG1-AS IN ( DUMMY ) 

DUMMY-NOSERAD/X1 
ANG2-ACOS (DUMMY) 

THETA3-PI-ANG1-ANG2 

DISTE-DISTC-NOSERAD+COS (THETA3 ) *NOSERAD 
DELANG-2 . 0 *THETA3 /( CRVGDY 1+ 1 ) 

DO  20  IY-1 , NOSGDY 


X ( IX+NOSGDX1 , IY , 1 ) -X ( IX+NOSGDX1-1 , I Y , I ) +DELDIST 
IF ( IY . LE. (INUMN/2) )THEN 

Y ( IX+NOSGDX1 , IY , 1) -CUMDELY 

Z (IX+NOSGDXl.IY, 1) “DISTD- (CUMDELY/TAN (THETA3) ) 

CALL  STRCH9 (DISTE, INUMN/2-I , IY, DELY, CIA, C1B, E1A, E1B, 
NOSGDX2 , IX) 

CUMDELY-CUMDELY+DELY 
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ELSEIF (IY . GE. (INUMN/2 ) +1 . AND. IY .LE. (INUMN/2)+CRVGDY141)THEN 

ICOUHT=ICOUNT+l 

DUMMY =THETA 3-1 COUNT  *  DELANG 

Y ( IX4NOSGDX1 , I Y ,  1 )  « ( DISTC-NOSERAD) +COS ( DUMMY ) *NOSERAD 
Z ( IX4NOSGDX1 , I Y , 1 ) «=SIN ( DUMMY) *NOSERAD 

ELSE 

Y ( IX+N0SGDX1 ,  IY,  1)  **Y  ( IX+N0SGDX1 , NOSGDY+l-I Y , 1 ) 

Z ( IX+NOSGDX1 , I Y , 1 ) =-Z ( IX+NOSGDX1 , NOSGDY+l-IY , 1 ) 

ENDIF 


20  CONTINUE 


PRINT  *, 'YES 3 ' 


BELOW  ARE  LISTED  SOME  MORE  VARIABLES  USED  IN  THE  PROGRAM: 

BODGDX1  -  THE  NUMBER  OF  SECOND  DELTA  WING  GRIDS  UP  TO  THE  RECTANGULAR 
SECTION  AND  NOT  INCLUDING  THE  FIRST  GRID  FROM  THE  NOSE 
NOSGDX  -  THE  TOTAL  NUMBER  OF  GRIDS  IN  THE  X  DIRECTION  OF  THE  FIRST 
DELTA  KING 

THETA4  -  THE  ANGLE  FORMED  BY  THE  SECOND  DELTA  KING 

DISTZU1  -  THE  DISTANCE  IN  THE  POSITIVE  Z  DIRECTION  FOR  THE  FIRST  NACA 
CROSS-SECTION  ENCOUNTERED  (DISTZU) 

DISTZD1  -  THE  SAME  DISTANCE  ON  THE  LOWER  SURFACE  (DISTZD) 

LEN4  -  THE  LENGTH  IN  THE  X  DIRECTION  FROM  THE  SECOND  DELTA  WING  TO 

THE  REAR  RECTANGULAR  SECTION 
LEN5  -  THE  TOTAL  LENGTH  OF  THE  SECOND  DELTA  WING 

THETA 5  -  THE  ANGLE  ON  THE  UPPER  SURFACE  FROM  THE  CENTERLINE  TO  THE 
FIRST  NACA  SECTION  ENCOUNTERED 
THETA6  -  THE  SAME  ANGLE  ON  THE  LOWER  SURFACE 

DISTA  -  THE  DISTANCE  IN  THE  Y  DIRECTION  WITH  A  NACA  CROSS-SECTION 

CURD  -  CHORD  LENGTH  OF  THE  NACA  SECTION 

XDIST  -  THE  POSITION  IN  THE  X  DIRECTION  ON  THE  NACA  AIRFOIL 


THE  SECTION  BELOW  GENERATES  THE  GRID  FOR  THE  SECOND  DELTA  WING 


NOSGDX=NOSGDXI+NOSGDX2 

CUMDIST=0.0 

DO  30  IX=NOSGDX+l , N0SGDX+B0DGDX1 
IXI=IX-N0SGDX1-N0SGDX2 

CALL  STRCH4 (LEN4 , BODGDX1 , IXI , DELDIST, DISTX2) 

CUMDIST=CUMDIST+DELDIST 

X(IX,1,I)=X(IX-1,1,1)+DELDIST 

Y(IX,1,1)«Y(IX-1,1,1) 

Z (IX, 1 , 1)«Z (IX-1 , 1 , 1) 

DISTE=LEN2+ ( CUMDIST/TAN (THETA4 ) ) -NOSERAD/SIN (THETA4 ) 
I COUNT =0 
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JCOTJNT=0 
KCOUNT=0 
CTJMDELV--C  .  0 

DO  30  IY«=2  ,  NOSGDY 

CALL  STRCH8  ( DISTE  ,  INUM/2 -1 ,  IY-1 , DELY , C2A ,  C2B ,  E2A,  E2 B , 
BCDGDX1 , IX-NOSGDX) 

CUMDELY-CUMDELY+DELY 
X(IX,IY,1)=X(IX-1,IY,1) +DELDIST 

IF (IY. LE . (INUM/2) ) THEN 

Y  (IX,  IY,  1)  *=Y (IX, IY-1, 1) +DELY 

ELSEIF (IY .GE. (INUM/2) +1 . AND. IY.LE. (INUM/2) + 
(CRVGDY-l)/2+l)THEK 

THETA9=ATAN ( ( Y ( IX , INUM/ 2 , 1 ) -Y ( IX , INUM/ 2- 1 , 1 ) ) / 

(Z (IX, INUM/2-1, 1) -2 (IX, INUM/2 , 1) ) ) 

DELANG=2 . 0*THETA9/ (CRVGDY+1) 

JCOUNT=JCOUNT+l 

YRAD=Z( IX, INUM/2, 1)/S IN (THETA9) 

Y  ( IX ,  I Y ,  1 )  <=Y  ( IX ,  INUM/2 , 1 )  +  (COS  (THETA9-DELANG* JCOUNT) 

*YRAD) -YRAD*C0S (THETA9) 


ELSE 

Y (IX, IY, 1) =Y ( IX , N0SGDY+1-IY , I) 

ENDIF 

CALL  NACA006 (DISTZUI , DISTZD1 , LEN5 , CUMDIST) 

THETA5=ATAN ( ( Z ( IX , 1 , 1 ) -DISTZU1 ) /LEN2 ) 

THETA6=ATAN ( ( Z ( IX , 1 , 1 ) +DISTZDI ) /LEN2 ) 

IF (IY. LT. ( (NOSGDY-1) /2 ) +1) THEN 

IF (Y (IX, IY , 1) . LE . LEN2 } THEN 

Z (IX, IY, 1)*Z (IX, 1, 1) -TAN (THETA5) *Y (IX, IY, 1) 

ELSEIF (Y  (IX,  IY,  1)  .GT. LEN2  .  AND.  IY. LE.  (INUM/2)  )THEN 

DISTA=Y (IX, IY , 1) -LEN2 
CHRD-LEN5-TAN (THETA4) *DISTA 
XDIST=CUMDIST-TAN (THETA4 ) *DISTA 

CALL  NACA006 (DISTZU, DISTZD, CKRD, XDIST) 

Z (IX, IY, 1) “DISTZU 

ELSE 

KCOUNT=KCOUNT  * 1 

Z (IX, IY, 1) “SIN (THETA9-DELANG*KCOUNT) *YRAD 
ENDIF 
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ELSEIF(IY.EQ. ( (NOSGDY-1) /2) +1)THEN 

KCOUNT=KCOUNT+ 1 
2 (IX, IY , 1) =0 . 0 

ELSE 

IF ( I Y . LE . INUM/2+CRVGDY+1) THEN 
KCOUNT=KCOUNT+ 1 

Z (IX, IY , 1) =SIN (THETA9-DELANG*KCOUNT) *YRAD 

ELSEIF (Y (IX,IY, 1) . GT . LEN 2 ) THEN 

DISTA=Y (IX, IY, 1) -LEN2 
CHRD=LEN5-TAN (THETA4 ) *DISTA 
XDIST=CUMDIST-TAN (THETA4 ) *DISTA 

CALL  NACA006 (DISTZU , DISTZD, CHRD, XDIST) 

Z ( IX, IY , 1) =DISTZD 

ELSE 

ICOUNT=ICOUNT+l 

Z(IX,IY,1)=Z (IX-1 ,NOSGDY , 1) +TAN (THETA6) *Y(IX,IY,1) 
ENDIF 
ENDIF 

30  CONTINUE 

PRINT  * , ' YES4 1 

THIS  SECTION  COMPUTES  THE  LAST  RECTANGULAR  SECTION  OF  THE  WING 


CUMDIST=0 . 0 
INUM2=NOSGDY-SIDGRD 
NOSERADl=NOSERAD*2 . 0 

DO  40  IX=NOSGDX+BODGDXl+l , N0SGDX+B0DGDX1+B0DGDX2 

IF(IX.GE.NOSGDX+BODGDX1+10)  INUM2-INUM2+3 
XI=( (LEN5-LEN4-NOSERAD) **2.0+ 

*  (Z (NOSGDX+BODGDX1 , IY , 1) -NOSERAD) **2.0)**0.5 
ANG1=ATAN ( ( Z (N0SGDX+B0DGDX1 , 1 Y ,  1) -NOSERAD) / 

*  (LEN5-LEN4-NOSERAD)  ) 

ANG2-ACOS (NOSERAD/X1 ) 

THETA7«=PI-ANG1-ANG2 
DELANG-THETA7/2 . 0 

DISTD=(LEN5-LEN4-NOSERAD+COS (THETA7) *NOSERAD) 

CUMDELY=0 . 0 

ICOUNT=0 

JCOUNT=0 

KCOUNT=0 

LCOUNT=0 
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MCOUNT=  0 

IXI=IX-N0SGDX1-N0SGDX2-B0DGDX1 

CALL  STRCH4 (DISTD, B0DGDX2 , IXI , DELDIST, DISTX3 ) 


CTJHDIST=CUMDIS?+D£LDIST 
DO  40  IY=1,N0SGDY 

XI= ( (LEN5-LEN4 -NOSERAD) **2 . 0+ 

*  (Z (NOSGDX+BODGDXl, IY,  1) -NOSERAD) **2 . 0) **0.5 
ANG1=ATAN ( (Z (NOSGDX+BODGDXl,  IY,  1) -NOSERAD) / 

*  (LEN5-LEN4-NOSEEAD)  ) 

ANG2=ACOS (NOSERAD/X1) 

THETA7=PI-ANG1-ANG2 
DELANG=THETA7/2 . 0 

DISTD= (LEN5-LEN4-NOSERAD+COS (THETA7) *NOSERAD) 

X(IX, IY , 1) =X (IX-1 , IY, 1) +DELDIST 
DISTE=LEN6-NOSERAD 

IF (IY.LE. ((NOSGDY-l)/2) +1) THEN 

IF (IY . LE . (INUM2/2) )TKEN 

Y  (IX,  IY,  1)  «=CUMDELY 

CALL  STRCH8 (DISTE, INUM2/2-1 ,  IY, DELY, C3A, C3B, E3A, E3B, 

*  BODGDX2 , IX-NOSGDX-BODGDX1 ) 

CUMDELY*=CUMDELY+DELY 

ELSEIF ( I Y . GT . INUM2 /2 . AND . 

*  IY.LE. INUM2/2+2) THEN 

icount=icount+i 

DELANG1«=(PI/4.u) 

Y(IX,IY,1)=Y(IX, INUM2/2 , 1) + 

*  COS ( PI/2. 0-ICOUNT*DELANG1)*NOSERAD 


ELSE 

Y(IX,IY,1)«Y(IX,IY-1,1) 

ENDIF 

ELSE 

J  COUNT=J  COUNT + 1 
Y (IX, IY , 1) =Y ( IX, IY-2*JCOUNT, 1) 

ENDIF 

IF (Y (IX, IY, 1) . LE.LEN2 . AND. IY. LE. INUM2/2 ) THEN 
Z(IX,IY,1)»Z(IX-1,IY,  1) -DELDIST/TAN (THETA7 ) 
IF ( Z ( IX , IY , 1 ) . LT.N0SERAD1)THEN 
Z(IX,IY,1)-N0SERAD1 
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ENDIF 


£LSEIF(Y (IX, IY( 1) . GT. LEN2 . AND. IY .  LE . INUM2/2 ) THEN 

DISTA=Y (IX, IY, 1) -LEN2 
CHRD=LEN5-TAN (THETA4 ) *DISTA 
XDIST=CHRD- (LEN5-LEN4 ) +CUMDIST 

CALL  NACA006 (DISTZU, DISTZD,  CHRD, XDIST) 

Z ( IX, IY , 1) =DISTZU 

IF(Z  (IX,  IY,  1)  .LT.KOSERADDTHEN 

Z (IX, IY , 1) =NOSERADl 

ENDIF 

ELSEIF (IY . GT. INUM2/2 .AND. IY. LE. INUM2/2+2 ) THEN 
LCOUNT=LCOUNT+l 

Z ( IX , IY , 1 ) =Z ( IX, INUM2/2 , 1 ) -NOSERAD+ 

*  SIN (PI/2-LCOUNT*DELANGl) *NOSERAD 

ELSEIF (IY.GT.INUM2/2+5. AND. IY.LT. (NOSGDY-1) /2+1) THEN 

DIFF= ( (NOSGDY-1) /2+1) - (INUM2/2+2) 

DEL=Z (IX, INUM2/2+2 , 1)/DIFF 
Z (IX, IY , 1) =Z ( IX , IY-1, 1) -DEL 

ELSEIF  (IY.  EQ .  (  (NOSGDY-1)  /  2 )  -f  1)  THEN 

2 (IX, IY , 1) =0 . 0 

ELSE 

MCOUNT=MCOUNT+l 

Z(IX,IY,1)=-1.0*Z(IX,IY-2*MCOUNT,1) 

ENDIF 

40  CONTINUE 

PRINT  * , • YES5 ' 


THIS  SECTION  DOES  THE  REPEATING  OF  THE  TRAILING  EDGE  GRID  TO 
FORM  THE  SECTION  OF  THE  SURFACE  GRID  THAT  FORMS  THE  WAKE 


TOT=NOSGDX+BODGDXl+BODGDX2 
DO  50  IX=1 , WAKGRD 

CALL  STRCH5 ( WAKLEN , WAKGRD, IX , DELX ,  DISTX4 ) 
DO  50  IY=1 , NOSGDY 
ITOT=TOT+IX 

X  ( ITOT,  IY ,  1 )  *=X  (ITOT-1 ,  IY ,  1)  +DELX 
Y (ITOT, IY, 1) =Y (ITOT-ix, IY, 1) 
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Z (ITOT , IY , 1) =Z ( ITOT-ix, IY,  1) 


50  CONTINUE 

PRINT  *,'YES6' 


C 

C  PRINT  THE  DATA  TO  A  FILE 
C 

C  DO  11  1=40,42 

C  DO  11  1=1 , NOSGDX+BODGDX1+BODGDX2+WAKEGRD 

C  DO  11  J=1 , NOSGDY 

C  PRINT  *,I,J,X(I,J,1) ,Y(I,J,1) ,Z(I,J,1) 

C  11  CONTINUE 


IZ=1 
IX=1 
IY=163 
REWIND  3 

WRITE  ( 3  )  I Y— IX-*- 1 ,  NOSGDY  ,  IZ 

WRITE (3)  ( (X(I,J, 1) , I=IX, IY) ,J=1, NOSGDY)  , 

*  ( (Y(I , J, 1) ,I=IX,IY) ,J=1, NOSGDY) , 

*  ( (Z(I,J, 1) , I=IX , IY) ,J=1, NOSGDY) 

CLOSE (UNIT=10) 

CLOSE (UNIT=11) 

STOP 

END 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

THIS  IS  THE  SUBROUTINE  NACA006  AND  IT  EVALUTATES  THE  Z  DIRECTION  VALUES 
THAT  ARE  ASSIGNED  TO  THAT  SPECIFIC  CROSS-SECTION 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


SUBROUTINE  NACA006 (ZPOS , ZNEG , CHORD , XVAL) 

DIMENSION  XPC(26) ,YPC(26) 

DATA  XPC/0. 0,0. 5, 0.75, 1.25, 2. 5, 5. 0,7. 5, 10. 0,15. 0,20. 0,25. 0,30.0, 

*  35.0,40.0,45.0,50.0,55.0,60.0,65.0,70.0,75.0,80,0,85,0, 

*  90.0,95.0,100.0/ 

DATA  YPC/0. 0,0. 494, 0.596, 0.754, 1.024, 1.405, 1.692, 1.928, 2. 298, 

*  2.572,2.772,2.907,2.981,2.995,2.919,2.775,2.575,2.331, 

*  2.050,1.740,1.412,1.072,0.737,0.423,0.157,0.0/ 

X=XVAL/ CHORD* 100.0 
DO  10  1=1,25 

IF (X. GT. XPC ( I ) .AND. X. LT.XPC (1+1 ) ) THEN 
IVAL=I 
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GO  TO  99 


ELSE 

CONTINUE 

ENDIF 

10  CONTINUE 

99  VAL1=X-XPC ( IVAL) 

VAL2=XPC(IVAL+1) -XPC(IVAL) 
VAL3=YPC ( IVAL+1 ) -YPC { IVAL) 

ZVAL= ( VAL1*VAL3/VAL2 ) +YPC ( IVAL) 
ZPOS=ZVAL*CHORD/100 . 0 
ZNEG=-1 . 0*ZPOS 

RETURN 

EKE 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  THIS  SUBROUTINE  STRETCHES  THE  GRID  IN  THE  X  DIRECTION 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  STRCKi (DIA, IDIV, IA, DELDIST) 

YDIST=DIA/IDIV 

IF ( IA . LE . IDIV/2 ) THEN 

IB=IA 

IC=1 


ELSE 

IB=IDIV-IA 

IC=-1 

ENDIF 

IF ( IA . GE . IDIV/ 2 . AND . IA . LE . IDIV/2+ 1 ) THEN 

DELDIST” ( (DIA/2.0) **2 . 0- ( (IDIV/2-1) *YDIST) **2 . 0) **0 . 5 
ELSE 

DIST1=( (DIA/2.0) **2.0- (IB*YDIST)**2. 0) **0.5 
DIST2=( (DIA/2.0) **2.0- ( (IB-IC) *YDIST) **2.0)  **0.5 
DELDIST=ABS (DIST1-DIST2) 

ENDIF 

RETURN 
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END 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  THIS  IS  THE  SUBROUTINE  FOR  STRETCHING  IN  THE  Y  DIRECTION  C 

c  c 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


SUBROUTINE  STRCH2 (RAD, JDIV, IA, DELY) 

YDIST-RAD/JDIV 
IF (IA. EQ. 1) THEN 

DELY= ( (RAD* *2 . 0) - ( (JDIV-IA) *YDIST) **2 . 0) **0.5 
ELSE 

DIST1« ( (RAD**2 . 0) - ( (JDIV-IA) *YDIST) **2 . 0) **0 . 5 
DIST2=( (RAD* *2. 0)-( ( JDIV-IA+1) *YDIST) **2 . 0) **0. 5 
DELY=ABS ( DIST1-DIST2 ) 

ENDIF 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C 

C  THIS  SUBROUTINE  IS  TOR  LINEAR  STRETCHING  IN  THE  X  DIRECTION 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


SUBROUTINE  STRCH4 (DIA, IDIV, IA, DELDIST, XVAR) 

TOT=  0 . 0 
SUBTOT-l.O 

DO  10  1=1 , IDIV/2-1 

SUBTOT=XVAR*  SUBTOT 
TOT=TOT+SUBTOT 

10  CONTINUE 

SUETOT= ( DIA/ 2 ) / (TOT+ 1 ) 

IF ( IA . LE . IDIV/ 2 ) THEN 

IB=IA 

ELSE 
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IB=1A-  ( IDIV/2 ) 

ENDIF 

IF ( IB . EQ . 1 ) THEN 
DELDIST=SUET07 
ELSE 

DO  20  1=1, IB-1 

SUBTOT=XVAR  *  SUBTOT 
20  CONTINUE 

DELDIST=SUETC7 

ENDIF 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  THIS  SUBROUTINE  IS  FOR  LINEAR  STRETCHING  IN  THE  WAKE  X  DIR  C 
C  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  STRCK5 (DIA, IDIV, IA, DELDIST, XVAR) 

TOT=0 . 0 
SUBTOT* 1.0 

DO  10  1*1, IDIV- 1 

SUBTOT=XVAR* SUBTOT 
TOT=TOT+SUBTCT 

10  CONTINUE 

SUBTOT=DIA/ (TOT-1) 

IF(IA.EQ. 1) THEN 

DELDIST=SUBT07 

ELSE 

DO  20  1=1 , IA-1 

SUBTOT=XVAR*SUETOT 
20  CONTINUE 
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non 


DELDIST=SUBTOT 


ENDIF 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

THIS  SUBROUTINE  IS  FOR  LINEAR  STRETCHING  IN  THE  V  DIRECTION  C 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  STRCH6 (DIA, IDIV, IA,  DELDIST, XVAR1 , XVAR2 , INUM, IY) 

TOT =0.0 
SUBTOT-1. 0 

XVAR=( ( IY-1 . 0) / (INUM-1 . 0) ) * (XVAR2-XVAR1) +XVAR1 

DO  10  1=1 , IDIV-1 

SUBTOT=XVAR ‘SUBTOT 
TOT=TOT+SUBTOT 

10  CONTINUE 

SUBTOT=DIA/ (TOT+1) 

IB=IDIV-IA • 1 

IF ( IB. EQ . 1) THEN 

DELDIST=SUBTOT 

ELSE 

DO  20  1=1, IB-1 

SUBTOT=XVAR*SUBTOT 
20  CONTINUE 

DELDIST=SUBTOT 

ENDIF 

RETURN 

END 

************************************************************************ 

************************************************************************ 

*  THIS  IS  A  SUBROUTINE  FOR  EXPONENTIAL  STRETCHING  IN  THE  Y-DIRECTION  * 


140 


SUBROUTINE  STRCH7  ( DIA ,  IDIV ,  I A ,  DELDIST ,  CONI ,  CON2  ,  EX1 ,  EX2  ,  KX ,  KA) 

CON= (KA-1 . 0) / (KX-1 . 0) * (CON2-CON1) +CON1 

EX= (KA-1 . 0) / (KX-1 . 0) * (EX2-EX1) +EX1 

DIA1=C0N*DIA**EX 

D£LY=DIA1/IDIV 


DISTY1=IA*DELY 
DISTY2= (IA-1 . 0) *DELY 
VAL1= (DISTY1/CON) ** (1-0/EX) 

VAL2= (DISTY2/CON) **(1.0/EX) 

DELDIST=VAL1-VAL2 

RETURN 

END 

********************************************************************** 
********************************************************************** 
*  THIS  IS  A  SUBROUTINE  FOR  SINEUSOI DALLY  STRETCHING  THE  Y-DIRECTION  * 
********************************************************************** 
********************************************************************** 


SUBROUTINE  STRCH8  (DIA,  IDIV ,  IA , DELDIST ,  Cl ,  C2  , El ,  E2  , KX , KA) 

CON= (KA-1. 0)/ (KX-1. 0) *(C2-C1)+C1 
EX= (KA-1. 0)/ (KX-1. 0)*(E2-E1)+E1 
RAD45=. 785398163 

DELY=RAD45/IDIV 
DISTY1=IA*DELY+RAD4  5 
DISTY2= (IA-1 . 0) *DELY+RAD45 
VALl’=CON*SIN  (DISTY1)  **EX 
VAL2=CON*SIN (DISTY2) **EX 
DIST=VAL1-VAL2 

DELDIST=DIST*DIA/ (1.0-SIN (RAD4 5)  ) 

RETURN 

END 


********************************************************************** 

a********************************************************************* 

*  THIS  IS  A  SUBROUTINE  FOR  SINEUSOI DALLY  STRETCHING  THE  Y-DIRECTION  * 


SUBROUTINE  STRCH9 (DIA , IDIV, IA , DELDIST,  Cl ,  C2 , El , E2 ,  KX ,  KA) 

CON= (KA-1. 0)/ (KX-1. 0) * (C2-C1)+C1 
EX= (KA-1. 0)/ (KX-1. 0) * (E2-E1 ) +E1 
PI=4 . 0*ATAN (1.0) 
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RAD4  5=PI/4 . 0 
RAD60=PI/3.0 
RAD90=PI/2.0 

RADSTP= (RAD90-RAD60) / (KX-1 ) 

DELY=( ( (RADSTP* (KA-1) ) +RAD60) -RAD45) / (IDIV) 

DISTY1* (IA-1) *DELY+RAD45 
DISTY2=(IA) *DELY+RAD4  5 
VAL1=C0N*SIN (DISTY1 ) **EX 
VAL2=C0N*SIN (DISTY2 ) **EX 
DIST=VAL2-VAL1 

DELDIST=DIST*DIA/ (SIN (RAD45+DELY*IDIV) -SIN (RAD45) ) 

RETURN 

END 
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APPENDIX  F. 


ADDITIONAL  SOURCE  CODE 
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THIS  PROGRAM  REMOVES  THE  FIRST  THREE  POINTS  OF  THE  HEMISPHERE,  BRINGS 
THE  APEX  TO  A  POINT,  RENUMBERS  THE  GRID  POINTS  AND  FINALLY  DOUBLE  THE 
THICKNESS  OF  THE  SURFACE  GRID  IN  THE  Z-DIRECTION 


DIMENSION  X (14 0,240, 60) ,Y (140, 240, 60) , Z (140 , 24 0, 60) , 

*  XX (14 0,24  0, 60) , YY ( 14  0, 240 , 60)  ,22(140,24  0,60) 

READ ( 3 )  II, JJ, KK 

READ ( 3 )  ( ( (X( J , K, L)  , J=1 , II) , K=1 , JJ) , L=1 , KK) , 

*  ( ( ( Y ( J , K, L) ,J=1,II),K=1,JJ) ,L=1,KK) , 

*  ( ( (Z(J,K,L) ,J=1,II) ,K=1,JJ) ,L=1,KK) 


DO  10  1=1,11-3 

DO  10  J=1,JJ 

IF (I . EQ. 1) THEN 

XX(1,J,1)=0.0 
YY ( 1 , J , 1) =0 . 0 
22(1, J,l) =0.0 

ELSE 

XX ( I , J , 1 ) =X ( 1+3 , J , 1 ) 
YY(I,J,1)=Y(I+3,J,1) 
22(I,J,1)=2(I+3,J,1)*2.0 

ENDIF 

10  CONTINUE 
REWIND  3 

WRITE (3)  XI-3 , JJ , 1 

WRITE (3)  ( ( XX ( J , K , 1 ) ,J=l,II-3) ,K=1,JJ) , 

*  ( (YY(J,K,1) ,J=l,II-3) ,K=1,JJ) , 

*  ( (Z2(J,K,1) ,J=l,II-3) ,K=1, JJ) 


STOP 

END 
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c 

C  THIS  PROGRAM  AGAIN  DELETES  THREE  GRID  POINTS  FROM  THE  NOSE  OF  THE 
C  SURFACE  GRID.  THIS  PARTICULAR  PROGRAM  THEN  CLOSES  THE  APEX  OF  THE 
C  SURFACE  GRID  TO  A  POINT  AND  FINALLY  RENUMBERS  THE  GRID  FOR  FIELD  GRID 
C  GENERATION 
C 


DIMENSION  X ( 150 , 240,60) , Y(150 , 24 0 , 60) , Z (150 , 240 , 60) , 

*  XX (150, 24 0,60) , YY (150,240,60)  ,22(150,240,60) 

READ ( 3 )  II , JJ , KK 

READ ( 3 )  (( (X(J,K,L)  , J=1 , II) , K=1 , JJ) ,L=1,KK) , 

*  (( (Y(J,K,L)  ,J=1,I1) , K=1 , JJ) , L=1 , KK) , 

*  ( ( (Z(J,K,L) , J=1 , II) , K=1 , JJ) , L=1 , KK) 


DO  10  1=1,11-3 

DO  10  J=1,JJ 

IF ( I . EQ . 1 ) THEN 

XX(1,J,1)=0,0 
YY (1 , J , 1) =0 . 0 
ZZ ( 1 , J , 1) =0 . 0 

ELSE 

XX ( I , J , 1 ) =X ( 1+3 , J , 1 ) 
YY(I,J,1)=Y(I+3,J,1) 
2Z(I,J,1)=Z(I+3,J,1) 

ENDIF 

10  CONTINUE 
REWIND  3 

WRITE ( 3 )  II-3 , JJ, 1 

WRITE ( 3 )  ( (XX (J,  K,  1) , J=l, I X — 3 ) ,K=1,JJ) , 

*  ((YY(J,K,1) ,J=l,II-3) ,K=1,JJ) , 

*  ( (ZZ(J,K,1) ,J=l,II-3) ,K=1,JJ) 


STOP 

END 


145 


noon 


THIS  PROGRAM  READS  THE  FIELD  GRID  AND  JUST  CHANGES  THE  Z-DIRECTION 
VALUE  BY  DOUBLING  THE  THICKNESS 


DIMENSION  X( 150, 240 ,60) , Y (150 , 240 , 60) , Z (150 , 240 , 60) , 

*  XX (150,240,60) , YY (150,240,60) , ZZ ( 150, 240 , 60) 

READ ( 3 )  II , JJ ,  KK 

READ (3)  ( ((X(J,K,L)  ,J=1,II) ,K=1,JJ) ,L=1,KK) , 

*  ( ( ( V ( J, K, L) ,J=1,II) ,K=1,JJ) ,L=1,KK) , 

*  (((Z(J,K,L) , J=l, II) , K=1 , JJ) ,L=1,KK) 


DO  10  1=1,11 

DO  10  J=1,JJ 

IF ( I . EQ. 1) THEN 

XX ( 1 , J , 1 ) =X ( I , J , 1 ) 
YY ( 1 , J , 1) =0 . 0 
ZZ(1,J,1)=0.0 

ELSE 


XX ( I , J , 1 ) =X ( I , J , 1 ) 

YY ( I , J , 1 ) =Y ( I , J , 1 ) 
ZZ(I,J,1)=Z(I,J,1)*2.0 

ENDIF 

10  CONTINUE 
REWIND  3 

WRITE ( 3 )  II,JJ,1 

WRITE ( 3 )  ( (XX (J ,  K,  1) ,J=1,II) ,K=1, JJ) , 

*  ( ( YY ( J , K, 1) , J=1 ,11) , K=1 , JJ) , 

*  ( (ZZ ( J , K, 1) ,J=1,II) , K=1 , JJ) 


STOP 


END 
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THIS  PROGRAM  ALSO  DELETES  THE  FIRST  THREE  POINTS  OF  THE  HEMISPHERE  ON 
THE  SURFACE  GRID.  BUT  THIS  PROGRAM  ALSO  DOUBLES  THE  THICKNESS  OF  THE 
SURFACE  GRID  IN  THE  Z-DIRECTION  ONLY 


DIMENSION  X(150, 240,60) , Y ( 150 , 240, 60) , Z (150 , 240 , 60) , 

*  XX (150, 240, 60) , YY (150,240,60) ,ZZ(150,240,60) 

READ ( 3 )  II , JJ ,  KK 

READ ( 3 )  ( ( (X(J,K,L) , J=1 , II) , K-l, JJ) , L=1 , KK) , 

*  ( ( (Y(J,K,L) ,J=1,II) , K=l, JJ) ,L=1,KK) , 

*  (((Z(J,K,L) , J=1 ,11) ,K=1, JJ) , L=1 , KK) 


DO  10  1=1,11-10 

DO  10  J=1,JJ 

I F ( I . EQ . 1 ) THEN 

XX(1,J,1)=X(I+10,J,1) 

YY ( 1 , J, 1) =0 . 0 
ZZ(1,J,1)=0.0 

ELSE 

XX(I,J,1)=X(I+10,J,1) 

YY (I , J , 1) =Y (1+10 , J , 1) 
ZZ(I,J,1)=Z(I+10,J,1)*2.0 

ENDIF 

10  CONTINUE 
REWIND  3 

WRITE ( 3 )  II-10,JJ,1 

WRITE ( 3 )  ( (XX(J ,K, 1) , J=1 , 11-10) ,K=1,JJ) , 

*  ( (YY ( J, K, 1) ,J=1,II-10) , X=1 , JJ) , 

*  ( (ZZ(J,K,1) ,  J=1 ,11-10) , K“1 , JJ) 


STOP 


END 


noonn 


THIS  PROGRAM  READS  THE  FINISHED  SURFACE  GRID  AND  DELETES  THE  FIRST 
THREE  POINTS  THAT  WERE  CONSTRUCTED  BY  THE  HEMISPHERE  AND  THEN 
RENUMBERS  THE  GRID  FOR  USE  BY  THE  FIELD  GRID  GENERATOR 


DIMENSION  X(120,240,60) ,Y (120 , 240, 60) ,2(120,240,60) , 

*  XX (120,240,60)  ,YY(1.20, 240,60)  ,  ZZ  (120,  240, 60) 

READ ( 3 )  II , JJ , KK 

READ (3)  (( (X(J,K,L) ,J=1,II) , K=l, JJ) , L=1 , KK) , 

*  ( ( (Y (J ,K,L) ,J=1,II) , K=1 , JJ) , L=1 ,KK) , 

*  ( ( (Z(J,K,L) ,J=1,II) , K=1 , JJ) ,L=1,KK) 


DO  10  1=1,11-10 

DO  10  J=1,JJ 

IF (I . EQ . 1) THEN 

XX(1,J,1)=X(I+10,J,1) 

YY (1 , J , 1) =0 . 0 
ZZ ( 1 , J , 1) =0 . 0 

ELSE 

XX (I, J, 1)=X(I+10,J,1) 

YY ( I , J , 1 ) =Y ( If 10 , J , 1 ) 
ZZ(I,J,1)=Z(I+10,J,1) 

ENDIF 

10  CONTINUE 
REWIND  3 

WRITE ( 3 )  II-10,JJ,1 

WRITE ( 3 )  ( (XX(J,K, 1) ,J=1,II-10) ,K=1,JJ) , 

*  ( (YY(J,K,1) ,J=1,II-10) ,K*1,JJ) , 

*  ( ( ZZ ( J , K, 1 ) , J=1 , 11-10) , K=1 , JJ) 


STOP 

END 
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c 

C  THIS  IS  A  PROGRAM  THAT  MANUALLY  ADJUSTS  THE  FIELD  GRID  NOSE  REGION 
C  FOR  THE  CYLINDRICAL  GRID.  IT  REDISTRIBUTES  THE  X  VALUE  ALONG  THE 
C  SINGULARITY  OF  THE  NOSE.  IT  ALSO  RESIGNS  THE  Y  AND  2  VALUES  FOR  THE 
C  FIELD  GRID. 

C 


DIMENSION  X( 117 ,240, 64) , Y ( 117 , 240 , 64 ) , Z (117 , 240, 64) , 

*  XX (14,240,64) , YY (14,240,64) ,ZZ(14,240,64) , 

*  XXX (130, 240, 64) , YYY (130 , 240 , 64 ) , ZZZ (130, 24 0 , 64 ) 

OPEN  (UNIT=1 0 ,  FILE= ' ddwnos-.  in 1  ,  STATUS* '  OLD 1 ) 

READ (10, *)  DIA, XVAR 
READ (24)  II , JJ , KK 

READ  (24)  (  (  (X  ( J ,  K,  L)  ,  J=1 ,  II )  ,  K=1 ,  JJ)  ,  L=1 ,  KK)  , 

*  ( ( (Y(J,K,L) ,J=1,II) ,K=1,JJ) ,L=1,KX) , 

*  ( ( (Z  ( J ,  K,  L)  ,J=1,II)  ,K=1,JJ)  ,L*=1,KK) 

DO  10  1=1, JJ 

DO  10  J=1 , KK 

XX ( 1 , 1 , J ) =X ( 1 , 1 ,  J ) 

YY ( 1 , I , J) =Y ( 1 , I , J) 

ZZ(1,I,J)=Z(1,I,J) 

10  CONTINUE 

DO  13  1=2,14 

CALL  STRCH (DIA, 13 , 1-1, DELDIST,XVAR) 

DO  13  J=1,JJ 
DO  13  K= 1 , KK 

IF ( I . EQ. 2 ) THEN 

XX (I , J, K) =DELDIST 
ELSE 


XX (I, J,K)=XX(I-1,J,K)+DELDIST 
ENDIF 

13  CONTINUE 

DO  11  1=2,14 

DO  11  J=1,JJ 

DO  11  K=1 , KK 

IF (K.EQ. 1) THEN 

YY (I , J , K) =0 . 0 
ZZ (I , J , K) =0 . 0 
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ELSE 


YY(I,J,K)=Y(1,J,K)-Y(1,J,1) 

ZZ(I,J,K)=Z(1,J,K)-Z(1,J,1) 

ENDIF 


11  CONTINUE 

DO  17  1=1,130 

DO  17  J=1,JJ 

DO  17  K=1 , KK 

IF (I . LE. 13 ) THEN 

XXX (I ,  J,  K) =XX ( 15-1 , J, K) 
YYY (I , J , K) =YY (I , J , K) 
ZZZ(I,J,K)=ZZ(I,J,K) 

ELSE 


XXX (I, J,K)=X(I-13, J,K) 
YYY(I,J,K)=Y(I-13,J,K) 

ZZ2 (I, J,K)=Z(I-13,J,K) 

ENDIF 

17  CONTINUE 
REMIND  34 

WRITE (34 )  130 ,JJ , KK 

WRITE (34)  ( ( (XXX ( J , K, L) ,J=1,130) , K=1 , JJ) , L=1 , KK) , 

*  ( ( (YYY (J, K, L) , J=l, 130) ,K=1,JJ) ,L=1,KK) , 

*  ( ( (ZZZ (J , K, L) , J=1 , 130)  ,K=1,JJ) ,L=1,KK) 


STOP 


END 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
C  THIS  IS  THE  SUBROUTINE  THAT  REDISTRIBUTES  THE  GRID  POINTS  USING  A 
C  LINEAR  TYPE  DISTRIBUTION. 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


SUBROUTINE  STRCH(DIA, IDIV, IA, DELDIST, XVAR) 

TOT=0 . 0 
SUBTOT=l . 0 

DO  10  1=1, IDIV- 1 


SUBT0T=XVAR*SUF707 

TOT-TOT+SUBTOT 


150 


u  u 


10  CONTINUE 


20 


SUBTOT=DIA/ (TOT+1) 

IF (IA. EQ. 1 ) THEN 

DELDIST=0 . O-SUBTOT 
ELSE 

DO  20  1=1 ( IA-1 

SUBTOT=XVAR*SUBTOT 

CONTINUE 

DELDIST=0 . O-SUBTOT 
ENDIF 
RETURN 


END 


o  o  o  o  o 


THIS  IS  A  PROGRAM  THE  READS  THE  FINISHED  SURFACE  GRID  AND  CAN 

BUILD  A  FILE  OF  ANY  PARTICULAR  YZ -PLANE  CROSS  SECTION  THAT  IS  DESIRED. 

THE  DESIRED  PLANE  IS  CHOSEN  BY  CHANGING  THE  "I"  VARIABLE. 


DIMENSION  X< 117, 240, 60) ( Y (117 , 24 0 , 60)  ,  Z  (117 , 240 , 60) , 

*  XX (1,24 0,60)  ,.YY(  1,24  0, 60)  , ZZ  (1, 240 , 60) 

READ (24 )  II , JJ, KK 

READ (24 )  ( ( (X ( J , K, L) ,J=1,II) ,K=1,JJ) ,L=1,KK) , 

*  ( ( (Y(J, K, L) ,J=1,II) ,K=1,JJ) ,L=1,KK) , 

*  ( ( (Z ( J , K, L) ,J=1,II) , K=1 , JJ) ,L-1,KK) 


DO  10  J=1,JJ 

DO  10  K=1 , KK 

XX (1,J,K)=X(I,J,K) 
YY(1,J,K)-Y(I, J,K) 
ZZ(1,J,K)=Z(I,J,K) 

10  CONTINUE 


REWIND  37 

WRITE (37)  1 , JJ , KK 

WRITE (37 )  ( (XX( 1 , K, L) ,K=1,JJ) ,L=1,KK) , 

*  ((YY(l,K,t),K«l,JJ),L*l,KK), 

*  ((ZZ(1,K,L) ,K=1,JJ) ,L=1,KK) 


STOP 


END 
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